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INTRODUCTION 


Utilization  of  advanced  composites  has  been  motivated  by 
their  high  strength  and  stiffness  properties,  as  well  as  by  the 
cost  advantages  associated  with  the  fabrication  processes  for 
composite  structures.  As  the  number  of  applications  increases 
cind  as  the  design  expertise  improves,  the  capability  to  assess 
their  response  to  various  loading  and  environmental  conditions 
becomes  increasingly  important.  Thus,  it  is  necessary  to  inves- 
tigate the  influence  of  such  factors  as  the  frequency-dependent 
damping  characteristics  of  the  constituents,  and  the  degradation 
of  material  properties  as  a function  of  moisture  and  temperature. 
The  effect  of  material  anisotropy  can  be  important  for  the  problem 
of  flexure,  where  the  shear  and  flexural  deformations  may  be  of  a 
comparable  order.  In  addition,  since  the  shear  properties  and 
the  damping  characteristics  are  primarily  matrix  dependent,  it 
appears  that  the  combined  effects  of  damping  and  shear  deformations 
will  play  an  important  role  in  the  dynamic  response  of  composite 
materials. 

The  dynamic  response  of  composite  materials  has  been  inves- 
tigated with  the  objective  of  assessing  the  influence  of  the  geom- 
etry and  the  material  variables  upon  the  vibrational  frequencies 
and  dispersion  characteristics  in  composite  laminates.  These 
studies  are  an  essential  first  step  in  assessing  the  response  of 
a structure  subjected  to  more  complex  transient  and  aerodynamic 
loads  which  are  frequently  experienced  in  aerospace  structures.  A 
direct  and  important  consequence  of  aerodynamic  loading  on  struc- 
tures is  flutter,  which  can  be  defined  as  the  dynamic  instability 
of  an  elastic  body  in  an  airstream.  A form  of  oscillatory  insta- 
bility which  differs  from  the  more  conventional  lifting  surface 
flutter  is  panel  flutter,  which  is  the  interaction  between  the  mis- 
sile or  aircraft  surface  skin  panel  motions  and  aerodynamic  forces 
in  the  supersonic  regime. 


The  effects  which  arise  from  the  viscoelastic  and  the  resulting 
damping  behavior  and  the  combined  effects  of  damping  and  shear 
deformation  in  composites  have  not  been  investigated  very  exten- 
sively, particularly  from  the  point  of  view  of  flutter.  Hence, 
the  problem  studied  here  is  that  of  panel  flutter  of  fiber- 
reinforced  epoxy  matrix  composite  laminates.  Emphasis  is  placed 
upon  incorporating  into  the  analysis  an  appropriate  constitutive 
relation  for  the  composite  laminate.  This  requires  consideration 
of  coupling  effects  between  extension  and  bending,  and  bending  and 
twisting.  It  also  requires  the  consideration  of  shear  deformation 
effects  as  a function  of  laminate  construction  and  stacking 
sequence  and  the  frequency-dependent  dynamic  viscoelastic  behavior 
of  composite  laminates.  In  addition,  it  is  known  that  edge 
restraints  of  the  panel  tend  to  limit  strongly  the  panel  flutter 
amplitude.  This  limitation  often  causes  the  modes  of  structural 
failure  to  be  those  peculiar  to  fatigue,  rather  than  explosive 
fracture  of  the  skin  surface.  Therefore,  it  is  also  necessary 
to  investigate  the  interaction  between  the  environmental  degrada- 
tion of  composite  laminate  properties  and  panel  flutter.  Thus, 
this  study  can  appropriately  be  termed  as  the  effects  of  environ- 
ment, and  damping  and  coupling  properties  of  composite  laminates 
on  panel  flutter. 

In  keeping  with  the  objectives  of  this  investigation,  the 
effort  consisted  of  the  following  tasks  (these  tasks  are  shown 
schematically  in  figure  1) : 

(i)  Derivation  of  f iber-rreinforced  epoxy  matrix  lamina/ 
laminate  effective  complex  moduli  by  utilizing  the  correspondence 
principle  for  vibrations  of  viscoelastic  composites.  A brief 
description  of  the  correspondence  principle  and  the  analytical 
derivations  are  given  in  Appendix  A,  while  a discussion  of  damping 
and  the  lamina/laminate  complex  properties  appears  in  the  section 
"Damping  Effects  in  Composite  Laminates." 
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(ii)  Determination  of  composite  laminate  shear  correction 
factors  by  matching  the  thickness-shear  mode  cut-off  frequencies 
(corresponding  to  infinitely  large  wave  lengths)  as  obtained  from 
a "refined"  shear  deformation  laminated  plate  theory  and  the  exact 
elasticity  solution.  An  approach  has  also  been  evolved  to  estimate 
the  transverse  shear  coupling  stiffness.  The  analytical  methodology 
is  described  in  Appendix  B and  the  coupling  and  shear-deformation 
effects  as  well  as  the  shear-correction  factors  as  obtained  from 
the  matching  technique  are  discussed  in  the  section  "Shear- 
Deformation  and  Coupling  Effects." 

(iii)  Preliminary  investigation  for  extending  the  constitutive 
relations  developed  in  (i)  and  (ii)  to  consider  the  effects  of 
moisture  and  temperature,  salient  features  of  which  are  outlined  in 
the  section  entitled  "Moisture  and  Temperature  Effects." 

(iv)  Analysis  of  composite  laminate  panel  flutter  by  utilizing 
the  laminate  constitutive  relations  developed  in  (i)  and  (ii)  and 
linear  piston  theory  aerodynamics.  Since  the  main  oojective  of 
this  study  is  to  evaluate  effects  of  environment,  attention  is 
restricted  to  response  of  laminates  in  cylindrical  bending.  The 
governing  differential  equations  are  derived  and  a method  for  their 
solution  is  outlined  in  Appendix  C.  Effects  of  damping  and  shear- 
deformation  on  critical  flutter  parameter  for  some  laminates  are 
evaluated  in  the  section  "Panel  Flutter  of  Viscoelastic  Composite 
Laminates . " 

Based  on  the  above-mentioned  analytical  developments,  predic- 
tive tools  (specifically,  computer  programs)  have  been  developed 
for  the  determination  of  (i)  lamina/laminate  effective  complex 
properties  from  constituent  complex  properties,  (ii)  shear  cor- 
rection factors  as  a function  of  laminate  construction  and  stack- 
ing sequence,  (iii)  critical  flutter  parameter  for  viscoelastic 
laminates  using  the  results  of  (i)  and  (ii),  and  (iv)  transient 
moisture  distribution  in  a laminate. 

Finally,  conclusions  are  drawn  from  the  results  thus  obtained 
and  areas  for  future  investigations  are  identified. 


LAMINATE  CONSTITUTIVE  RELATIONS 


DAMPING  EFFECTS  IN  COMPOSITE  LAMINATES 

A vibrating  structure,  such  as  a composite  laminate  panel, 
has  at  any  instant  some  kinetic  energy  and  some  strain  (or  potential) 
energy.  The  kinetic  energy  storage  is  associated  with  the  mass 
and  velocity  at  that  instant,  and  the  strain  energy  storage  with 
the  stiffness.  In  addition,  the  structure  also  dissipates  some 
energy  as  it  deforms.  It  is  this  energy  dissipation  - or,  more 
precisely,  this  conversion  of  mechanical  energy  into  thermal  energy  - 
that  is  called  damping.  Unlike  mass  and  stiffness,  damping  does 
not  refer  to  a unique  physical  phenomenon,  and  hence,  damping  is 
much  more  difficult  to  define.  There  are  a multitude  of  damping 
mechanisms  such  as;  interface  friction,  fluid  viscosity,  turbulence, 
acoustic  radiation,  eddy  currents,  magnetic  hysteresis,  and  mech- 
anical hysteresis  (also  called  internal  friction  or  material  damp- 
ing) . In  order  to  analyze  or  predict  the  damping  of  a given  struc- 
ture, one  should  ideally  take  into  account  all  possible  damping 
mechanisms;  however,  in  practical  cases,  one  or  two  mechanisms 
generally  predominate. 

The  primary  effects  of  damping  are  (i)  reduction  of  vibration 
amplitudes  at  resonances,  with  attendant  reductions  in  stresses, 
structural  fatigue,  and  sound  radiation;  (ii)  rapid  decay  of  free 
vibrations;  (iii)  attenuation  of  structureborne  waves  propagating 
along  the  panel,  i.e.  reduced  conduction  of  vibratory  energy  along 
the  panel;  and  (iv)  increased  sound  isolation  (transmission  loss). 

All  these  effects  are  generally  beneficial  from  the  standpoint  of 
noise  and  vibration  control.  Since  the  problem  of  airflow  over  a 
thin  panel  involves  both  the  noise  and  vibration  aspects,  damping 
is  expected  to  have  an  important  effect  on  panel  flutter. 

Lamina/Laminate  Complex  Properties 

The  reason  for  the  occurrence  of  damping  in  fiber-reinforced 
composites  is  the  dissipative  nature  of  epoxy  matrix.  The  fiber 
and  matrix  may  be  characterized  as  linear  viscoelastic  materials. 


The  problem  of  damping  in  composite  materials  has  been  investigated 
both  analytically  (refs.  1-3)  and  experimentally  (refs.  4-6).  If 
the  fiber  and  matrix  complex  moduli  (storage  and  loss  moduli)  are 
Jcnown  from  experiments  as  a function  of  frequency,  it  is  a relatively 
straightforward  task  to  obtain  the  lamina/laminate  complex  moduli 
by  utilizing  the  dynamic  viscoelastic  correspondence  principle. 

This  principle  states  that  "the  effective  complex  moduli  of  a visco- 
elastic composite  are  obtained  by  replacing  the  phase  elastic  moduli 
by  the  phase  complex  moduli  in  the  expressions  for  the  effective 
elastic  moduli  of  identical  phase  geometry."  If  the  effective 
elastic  moduli  are  of  a form  such  that  they  are  explicitly  obtained 
from  the  constituent  properties,  as  in  references  7 and  8,  the 
dynamic  correspondence  principle  is  particulary  simple  to  apply. 

The  logic  for  the  correspondence  principle  for  vibrations  of 
viscoelastic  composites  is  reconstructed  from  reference  1 in 
Appendix  A.  The  expressions  for  the  lamina  complex  moduli  as  a 
function  of  phase  complex  moduli  (from  the  elastic  composite  cylin- 
der assemblage  model  of  references  7 and  8)  and  the  laminate  (see 
fig.  2)  complex  stiffnesses  are  also  given  in  the  Appendix.  Results 
for  the  lamina/laminate  complex  properties  obtained  from  these 
expressions  are  discussed  subsequently. 

Discussion  of  Results 

The  constituent  properties  used  to  obtain  the  results  given 
in  this  section  are  tabulated  in  table  1.  Note  that  the  fiber 
properties  are  assumed  to  be  elastic  and  frequency  independent. 

The  frequency-dependent  Young's  moduli  of  epoxy  are  obtained  from 
reference  2 and  extrapolated  for  frequencies  greater  than  10  . 

The  loss  factor  is  seen  to  increase  with  the  frequency.  In  general, 
however,  the  loss  factor  frequency  curve  passes  through  peak(s) 
and  trough (s)  at  certain  frequencies.  A change  in  temperature  is 
likely  to  alter  the  location  of  these  peaks  and  peak  values.  A 
good  amount  of  data  on  damping  characteristics  of  metals  is  avail- 
able in  literature.  For  epoxy  matrices,  however,  such  datei  is 
scarce.  Poisson's  ratios  for  various  frequencies  given  in  table  1 
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are  calculated  on  the  assumption  that  the  matrix  material  is 
elastic  in  dilation  and  the  bulk  modulus  is  .3824  x 10^  psi,  a 
constant  for  all  frequencies. 

The  variation  of  Graphite/ epoxy  lamina  complex  Young's  and 

shear  moduli  and  Poisson's  ratio  as  a function  of  frequency  is 

illustrated  in  figures  3,  4,  and  5,  respectively.  In  figure  3, 

the  real  part  of  the  lamina  axial  modulus  is  insensitive  to  the 

frequency.  This  is  to  be  expected  because  the  fiber  properties 

are  assumed  to  be  frequency  independent.  The  transverse  lamina 

properties  (both  real  and  imaginary),  the  imaginary  part  of  the 

axial  modulus,  and  the  shear  moduli  in  figure  4,  however,  follow 

the  matrix  property  variation  with  frequency.  The  complex  laminate 

stiffnesses  ”^11'  *^44  45  ^ function  of  frequency  are  shown 

in  figures  6 through  8 for  [0/90]  , [0^/±45]  , and  [0/±45/90] 

s z s ^ s 

Graphite/epoxy  laminates.  For  the  complex  extensional  and 

bending  stiffnesses,  only  the  imaginary  parts  are  functions 

of  frequency  while  for  the  complex  transverse  shear  stiffnesses 
(A44  , the  variation  of  both  the  real  and  imaginary  parts  is 

governed  by  the  matrix. 

In  summary,  the  complex  stiffnesses  for  a composite  lamina/ 
laminate  are  obtainable  in  a simple  manner  from  the  phase  complex 
moduli.  However,  at  the  present  time,  data  are  lacking  for  the 
properties  of  the  fiber  and  matrix  at  higher  values  of  frequencies. 

SHEAR-DEFORMATION  AND  COUPLING  EFFECTS 

A number  of  theories  have  been  developed  in  the  literature  to 
characterize  the  mechanical  response  of  fiber-reinforced  composite 
laminates.  These  range  from  the  classical  laminated  plate  theory 
which  employs  the  Kirchhof f-Love  hypothesis  to  the  solutions  of  the 
exact  elasticity  equations.  Since  in  composites  the  ratio  of  the 
Young's  modulus  to  the  shear  modulus  can  be  very  large,  transverse 
shear  deformations  in  plates  of  such  materials  can  be  quite  signifi- 
cant. For  this  reason,  laminated  plate  theories  based  on  Kirchhoff- 
Love  hypotheses  (refs.  9-10)  are  not  always  adequate  for  predicting 


structural  response.  Failure  to  consider  shear  deformation  may 
lead  either  to  overestimates  of  the  frequencies  of  vibration  or 
underestimates  of  static  deflections  for  the  smaller  plate  dimension 
to  the  plate  thickness  ratios  of  less  than  or  equal  to  10.  Shear 
deformations  have  been  considered  in  plate  theories  (refs.  11  '.6) 
in  a manner  similar  to  that  in  Timoshenko- type  beam  theory  (ref.  17) 
by  introducing  shear  correction  factors.  The  laminated  plate 
response  can  also  be  predicted  by  higher  order  theories  such  as 
the  elasticity  theory  or  the  "continuum"  theory.  Solutions  to 
the  static  and  dynamic  problems  have  been  obtained  (refs.  18-20) 
by  utilizing  the  "exact"  elasticity  theory.  The  scope  of  this 
analysis  is,  however,  limited  because  of  the  lengthy  and  tedious 
computations  involved.  In  the  "continuum  theory"  (refs.  21-22) 


for  a layered  medium  consisting  of  alternating  plane  parallel  layers 
of  homogeneous  isotropic  materials,  the  local  layer  rotations  are 
assumed  distinct  and  the  transition  from  the  pair  of  alternating 
layers  to  the  continuum  is  achieved  by  a smoothing  operation.  This 
theory  has  also  been  referred  to  as  the  "effective  stiffness" 
theory.  The  drawback  of  the  "effective  stiffness"  theory  is  that 
its  applicability  to  "real-world"  composite  laminates  is  limited. 

The  coupling  effects  in  composite  laminates  arise  out  of  the 
extension-bending  and  bending-twisting  stiffnesses.  For  example, 
the  stiffness  B^j(i,j  = 1,2,6)  in  equation  (B17)  in  Appendix  B 
represents  the  coupling  phenomenon  between  stretching  and  bending. 


For  a symmetrical  layup  of  the  laminate, 
eliminated.  For  unsymmetrical  laminates 


B.  . 
ID 

with 


0 and  coupling  is 
alternating  +6®/-0® 


layers  of  equal  thickness,  it  can  be  shown  that 


' ®26)  = 


' ^26^ 


(1) 


Thus  by  increasing  the  number  of  layers  m in  a laminate  (fig.  2 ) , 
B^g  and  become  smaller  in  magnitude  and  the  bending  extension 

coupling  is  reduced.  This  is  because  the  laminate  becomes  more 
homogeneous  with  increasing  m.  The  effect  of  coupling  is  to  re- 
duce the  natural  frequency  of  a laminate  and  the  degree  of  the 
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effect  depends  upon  the  number  and  orientation  of  the  layers  in 
the  laminate.  Implications  of  the  coupling  effect  will  be  evident 
on  panel  flutter  by  virtue  of  the  change  in  the  natural  frequency. 
Also,  in  the  presence  of  uniform  in-plane  loads  (prestressed  panels) , 
the  coupling  effect  will  assume  particular  significance. 

Shear  Correction  Factors  for  Composite  Laminates 

As  mentioned  earlier,  the  need  for  the  consideration  of  shear 
deformation  in  composite  laminates  arises  from  the  fact  that  the 
material  is  "soft"  in  shear.  Also,  because  the  modes  of  interest 
for  the  panel  flutter  problem  are  primarily  the  flexural  and  the 
thickness  shear  modes,  it  is  logical  to  consider  a theory  which 
would  include  only  the  pertinent  deformations  and  eliminate  the 
unwanted  modes  or  deformations  associated  with  it  such  as  the  ex- 
tensional  modes.  Besides,  it  has  also  been  demonstrated  in  refer- 
ence 20  that  the  shear  deformation  theory  shows  good  correlation 
with  the  exact  elasticity-solution  for  the  first  few  modes  of 
vibration.  Associated  with  this  theory  is  the  shear  correction 
factor,  k,  occurring  with  the  transverse  shear  stiffnesses 
Aij(i,j  = 4,5).  These  correction  factors  are  obtained  by  matching 
the  results  of  the  approximate  plate  theory  with  those  obtained 
from  more  accurate  solutions,  either  for  the  static  (energy  formu- 
lation; ref. 11  ) or  the  dynamic  problem  (thickness-shear  mode 
frequency  matching  technique;  ref.  12).  For  the  static  case  of 
isotropic  homogeneous  plates  (ref.  11),  this  shear  correction  fac- 
tor is  obtained  to  be  k = 5/6,  while  for  the  dynamic  case  (ref.  12), 
2 

it  is  k = TT  /12.  It  is  also  well  known  that  the  shear  correction 

2 

factor  for  an  orthotropic  homogeneous  plate  is  the  same  (=7t  /12)  as 
that  for  an  isotropic  plate.  The  approach  generally  followed  for 
composite  laminates  is  to  use  a constant  value  of  the  shear  coeffi- 
cient of  the  shear  correction  factor  for  different  laminates  ir- 
respective of  the  laminate  construction  and  stacking  sequence. 

There  is,  however,  considerable  uncertainty  associated  with  using  a 
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constant  value  of  k for  all  transverse  shear  stiffnesses  (A..,  A._ 

44  4 5 

and  A^^)  for  all  laminates  and  stacking  sequences.  It  was  noted 

in  reference  20  that  k should  be  different  for  the  deformations 

associated  with  y and  y ^nd  that  it  should  be  also  a function 

' yz  xz 

of  fiber  orientation  and  stacking  sequence. 

Referring  again  to  equation  (B17)  in  Appendix  B,  the  extension- 
bending and  bending  (D^j)  stiffnesses  are  functions  of  stacking 

sequence  while  the  transverse  shear  stiffnesses  (A^ ^ ; i , j=4 , 5)  are  not. 
Thus,  a constitutive  theory  should  be  developed  which  is  accurate 
enough  to  consider  all  the  stacking  sequence  effects  but  simple 
enough  to  be  used  in  a form  suitable  for  structural  dynamics.  Hence, 
an  approach  has  to  be  formulated  by  which  the  transverse  shear  stiff- 
ness coefficients  which  are  ordinarily  independent  of  the  stacking 
sequence  can  be  obtained  as  a function  of  stacking  sequence.  The 
only  manner  in  which  this  can  be  achieved  by  utilizing  a simple 
shear  deformation  theory  is  to  develop  a procedure  to  obtain  the 
shear  correction  factors  associated  with  these  stiffnesses  as  a 
function  of  stacking  sequence  and  laminate  construction.  While  this 
problem  has  been  considered  in  references  3,  15,  and  16,  a system- 
atic approach  to  study  the  effects  of  laminate  construction  and 
stacking  sequence  is  lacking. 

In  this  study,  the  shear  correction  factors  are  obtained  by 
utilizing  the  exact  elasticity  theory  and  considering  only  the 
thickness-shear  modes  of  vibration.  Cut-off  frequencies  are  obtained 
for  infinitely  long  wave  lengths  in  the  limit  and  these  frequencies 
are  compared  with  the  cut-off  frequencies  for  the  same  thickness 
shear  mode  of  vibrations  obtained  from  the  shear  deformation  theory. 
Since  the  cut-off  frequencies  will  theoretically  not  match  because 
of  the  approximate  nature  of  the  shear  deformation  theory,  they  are 
forced  to  match  by  the  inclusion  of  the  shear  correction  factor  in 
the  shear  deformation  theory.  Since  the  exact  elasticity  solution 
accounts  not  only  for  the  laminate  construction,  but  also  for  the 
stacking  sequence,  it  is  evident  that  the  shear  correction  factors 
so  obtained  will  indeed  be  functions  of  laminate  construction  and 


stacking  sequence.  In  this  regard,  two  approximate  shear  deformation 
theories  have  been  considered;  namely,  the  one  which  does  not  account 
for  the  matching  of  interlaminar  stresses  at  the  interfaces  by  con- 
sidering a constant  transverse  shear  strain  through  the  thickness  of 
the  laminate  (refs.  13  and  14),  and  the  one  which  considers  piece- 
wise  continuous  transverse  shear  strains  through  the  laminate  in 
an  average  fashion  by  matching  the  interlaminar  stresses  at  the 
interface  (ref.  23) . It  is  observed  that  the  shear  correction  fac- 
tors as  obtained  from  each  of  these  approximate  shear  deformation 
theories  differ  significantly,  though  the  products  of  the  shear 
corrections  factors  and  the  corresponding  transverse  shear  stiff- 
nesses (defined  by  the  matrix  i,j=4,5)  are  identical.  It  is 

interesting  to  note,  however,  that  the  shear  correction  factors 
obtained  from  the  latter  theory  for  a large  number  of  layers 
approach  the  value  for  a homogeneous  material,  while  those  from 
the  former  theory  do  not. 

The  transverse  shear  coupling  stiffness  as  obtained  in 

equation  (B17)  in  Appendix  B,  is  zero  for  orthotropic,  symmetric,  and 

[±6]  laminates  (±6  laminates  having  equal  percentages  of  +0  and  -0 

layers).  Thus  the  coupling  stiffness  ^ ^^45^45  ^45^45^  also 

vanishes  for  these  laminates.  In  reality,  however,  the  transverse 

shear  force  Q will  induce  both  y and  y shear  strains,  although 
X xz  yz 

the  magnitude  of  the  latter  will  be  smaller  than  that  of  the  former. 
The  actual  coupling  is  expected  to  reduce  for  laminates  which  are 
more  homogeneous.  In  the  present  investigation,  an  approach  has 
been  developed  to  assess  the  magnitude  of  this  coupling. 

The  analytical/numerical  procedure  to  obtain  the  shear  correc- 
tion factors  by  matching  the  thickness-shear  mode  cut-off  frequen- 
cies obtained  from  the  exact  elasticity  solution  and  the 
shear-deformation  theory  is  described  in  Appendix  B.  In  addition, 
the  method  for  the  evaluation  of  the  transverse  shear  coupling  term 
has  been  outlined  in  the  Appendix. 
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Discussion  of  Results 


The  shear  correction  factors  obtained  by  the  methods  described 

above  and  in  Appendix  B are  shown  in  figures  9 through  24  for  various 

mid-plane  symmetric  laminates  with  the  fiber  direction  of  the  0° 

layer  coinciding  with  the  x direction  (see  fig.  2) . The  results 

are  reported  for  various  values  of  the  ratio  r G23/Gj_2(r<  D 

where  Gj^2  *^23  ^~^13  axial  and  transverse  shear  modulus, 

respectively,  of  each  layer.  In  figure  9,  the  shear  correction 

factors  for  a [0/90]  laminate  are  based  on  the  shear-deformation 

s 

theory  of  reference  14  which  assumes  a constant  transverse  shear 
strain  across  the  thickness.  It  is  evident  from  the  figure  that 
for  increasing  number  of  layers,  the  correction  factors  and  k^^ 

(k^^  = 0 for  this  laminate)  do  not  approach  the  value  tt‘'/12  for  a 
homogeneous  plate.  This  phenomenon  has  also  been  observed  in 
reference  16  for  the  static  solution.  If  a more  "refined"  shear 
deformation  theory  is  used  which  considers  a constant  value  of  inter- 
laminar shear  stress  through  the  thickness  (ref.  23) , it  is  seen 

I I 

that  the  shear  correction  factors  k^^  and  k^^  for  the  [0/90] ^ 
laminate  (fig.  10)  indeed  approach  the  value  for  the  homogeneous 
plate.  As  noted  in  Appendix  B,  it  is  of  no  consequence  as  to  which 
shear  correction  factors  are  used  (kj^j  ^ij^'  provided  that  the 
shear  stiffness  or  kj^^A^^)  are  identical. 

Results  for  mid-plane  symmetric  angle-ply  laminates  (±0, 

0 = 15°,  30°,  40°,  and  45°)  are  given  in  figures  11-18  and  are 

I 

summarized  in  figures  19  and  20.  Since  A^^  is  equal  to  zero  for 
symmetric  laminates,  the  correction  factor  k.^.  cannot  be  calculated 

4 !>  I, 


and  does  not  have  any  significance. 


A quantity  k^^  is  defined  with 


the  help  of  which  the  coupling  stiffness  K.j-  can  be  calculated  as 

" 1/2  ^ ^ 

(K^^K^g)  ' . These  results  show  that  the  coupling  shear 


K 


45 


= k 


45 


correction  factors  k^^  approach  zero  and  k^^  and  k^^  approach 
7r^/12  with  increasing  m,  as  in  the  case  of  cross-ply  laminates. 

It  should  be  pointed  out  that  results  for  r > 1 can  be  easily 

t • 

obtained  from  the  results  given  here  by  interchanging  k^^  and  k^^ 
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r 


II 

for  the  case  G_-/G, _ = 1/r  and  leaving  unchanged.  For  r = 1, 

I I 


44 


= kcc  = IT  /12  and  k._  = 0, 
5d  45 


For  m > 8 , the  correction  factors 
2 


^44'  ^55  differ  appreciably  from  it^/12  for  all  the  laminates 

considered  provided  that  r is  not  too  small.  Hence,  the  value 

2 ' 
of  IT  /12  can  be  used  in  the  plate  theory  if  are  obtained  in 

a manner  as  outlined  in  Appendix  B and  reference  23. 


Correction  factors  k^^,  k^^,  and  k^g  for  a [0/+45/-45/90] ^ 

laminate  with  equal  percentages  of  0°  and  90®  layers  are  shown  in 

figure  21  for  different  percentages  of  ±45°  layers  and  various 

values  of  the  ratio  r.  Similar  plots  for  [0/±45]  , [±45/0/90]  , 

s s 

and  (±45/0)  laminates  are  given  in  figures  22-24.  Figures  21-24 

® M 

show  that  the  coupling  term  k^^  is  extremely  small  when  the 
volume  of  ±45°  layers  (equal  volume  of  +45  and  -45)  is  less  than 
0.5.  In  general,  for  laminates  with  a large  number  of  0°,  90°, 

I 

and  ±0  layers  arranged  in  a periodic  fashion,  the  term  k^^  is 
likely  to  be  small. 

) I 

Correction  factors  k^^ , k^^  for  antisymmetric  laminates  con- 
structed with  a sequential  arrangement  of  +0  and  -0  layers  are 
plotted  against  0 for  various  values  of  r in  figure  25.  It  may 
be  noted  that  for  r =.2,  a discontinuity  exists  in  the  plots  for 
a value  of  0 = 0^,  where  0q  depends  on  the  total  number  of  layers 
m,  but  is  very  close  to  12°.  Near  the  discontinuity,  the  cut-off 
frequencies  corresponding  to  the  first  u-dominated  mode  (anti- 
symmetric) and  the  second  v-dominated  mode  (symmetric)  are  close 
to  each  other.  Also,  they  do  not  coalesce  at  0 = 0^.  For  6 < 0q, 

the  larger  frequency  corresponds  to  the  u-dominated  mode  and  the 
lower  one  is  for  the  v-dominated  mode.  As  0 passes  through  0^, 
the  roles  are  reversed.  Near  ® the  antisymmetric 

u-displacement  field  is  strongly  coupled  with  symmetric  v displace- 
ments. The  laminated  plate  theory  cannot  predict  any  symmetric 

I 

displacement  fields.  The  correction  factors  k^^  are  calculated 
so  as  to  predict  the  frequency  corresponding  to  the  u-dominated 
mode,  but  the  usefulness  of  such  factors  near  0 = 0^  is  limited. 


i 
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The  discontinuity  is  pronounced  when  the  number  of  layers  m is 
small  and  cannot  be  identified  for  m = 8.  Fortunately,  for 
commonly  used  laminates,  m is  large  and  r is  of  the  order  of  0.5 
or  higher,  and  the  laminated  plate  theory  can  be  used  to  predict 
the  plate  response  with  reasonable  accuracy  in  such  cases. 

However,  when  G22  sind  G^2  of  different  order,  the  laminated 

plate  theory  should  be  used  with  caution.  The  coupling  between 
antisymmetric  u and  v displacements  is  extremely  small  in  these 
laminates,  and  the  coupling  stiffness  can  be  taken  equal  to  zero. 

It  should  be  pointed  out  that  as  the  number  of  layers  is 

' ' 2 

increased,  quickly  approach  fr  /12. 

Laminated  plates  can  be  constructed  with  various  stacking 
sequences  of  layers  of  different  orientations.  Correction  factors 
for  only  a limited  number  of  such  laminates  are  presented  here. 

The  results  demonstrate  the  usefulness  of  the  approach  and  indicate 
some  general  trends  which  will  be  of  value  in  design  of  such  plates. 
In  this  regard,  it  should  be  mentioned  that  for  the  application 
considered  in  this  study  (panel  flutter  of  viscoelastic  laminates) , 
shear  moduli  of  each  layer  are  frequency-dependent  complex  quan- 
tities. To  determine  the  effects  of  shear  deformation  of  such 
laminates,  one  may  consider  the  problem  of  wave  propagation  in  a 
heterogeneous  viscoelastic  medium.  Alternatively,  one  can  take 
recourse  to  the  dynamic  correspondence  principle  (ref.  1)  provided 
analytical  expressions  for  correction  factors  are  available.  The 
solution  presented  here  does  not  yield  closed-form  analytical 
expressions.  However,  for  a given  stacking  sequence  of  layers, 
the  correction  factors  depend  only  on  the  ratio  r ^22^^12' 
Therefore,  if  the  imaginary  part  of  this  ratio  in  a viscoelastic 
case  is  small  and  the  real  part  is  fairly  constant  over  a frequency 
range  near  the  cut-off  frequencies,  the  values  given  here  can  be 

used  as  reasonable  estimates  in  the  viscoelastic  case  as  well. 

2 

In  particular,  the  limiting  value  of  tt  /12  for  a laminate  made  of 
a large  number  of  layers  arranged  in  a periodic  fashion  holds  both 
for  the  elastic  and  viscoelastic  cases. 


MOlS>rbRE  AND  TEMPERATURE  EFFECTS 

An  important  aspect  associated  with  the  laminate  constitutive 
equations  is  the  effect  of  environment.  Since  the  dynamic  response 
characteristics  of  the  composite  panel  must  be  known  for  its  entire 
lifetime,  the  effects  of  changing  environment  and  material  wearout 
or  degradation  should  be  considered.  Recent  developments  in  the 
application  of  polymeric  composites  to  aircraft  structures  have  re- 
vealed significant  changes  in  the  mechanical  properties  and  dimen- 
sions of  these  materials  due  to  the  effects  of  combined  absorbed 
moisture  and  thermal  environment  (refs.  24-29).  Material  degrada- 
tion will  primarily  manifest  itself  in  the  form  of  property  (stiff- 
ness as  well  as  strength)  losses.  It  has  been  observed  that  absorbed 
moisture  coupled  with  thermal  spikes  adversely  influences  the  matrix 
controlled  properties.  Additionally,  absorbed  water  molecules  may 
act  as  a plasticizer  and,  therefore,  reduce  the  glass  transition 
temperature  of  the  matrix,  i.e.,  the  temperature  at  which  there  is  a 
strong  loss  of  matrix  modulus  and  with  it,  a corresponding  loss  of 
composite  strength  and  a significant  increase  in  material  deimping. 

For  the  panel  flutter  problem,  the  loss  in  the  transverse 
shear  and  flexural  stiffnesses  due  to  interlaminar  degradation 
will  affect  critical  flutter  speeds.  Also,  the  alteration  in 
the  glass  transition  temperature  due  to  humidity  and  temperature 
will  alter  the  damping  characteristics  of  the  matrix  (refs.  29-31) 
and  hence , the  laminate . 

The  analysis  approach  for  the  determination  of  the  effects  of 
moisture  and  temperature  on  the  laminate  constitutive  equations  is  to 
determine  the  launina  dif fusivities  from  constituent  dif fusivities 
using  the  conductivity  analogy  and  composite  cylinder  assemblage 
results.  Having  obtained  the  lamina  dif fusivities  (and  these  dif- 
fusivities  could  be  functions  of  moisture  content  and  temperature) , 
it  is  possible  to  calculate  the  laminate  through-the-thickness 
transient  or  steady  state  moisture  distribution  in  the  leuninate  from 
the  one-dimensional  diffusion  equation  (refs.  32-33).  If  the  lamina 
dif fusivities  are  independent  of  the  moisture  content  and  if  the 
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problem  is  solved  for  a given  temperature,  which  is  constant  through- 
out the  laminate,  the  solution  of  the  diffusion  equation  is  rather 
straightforward.  However,  if  one  considers  the  dif f usivities  as  a 
function  of  moisture  and/or  temperature  and  if  there  is  a tempera- 
ture gradient  across  the  laminate  then  the  solution  has  to  be  ob- 
tained by  numerical  methods.  The  reason  for  obtaining  the  transient 
state  moisture  content  in  a laminate  is  that  a laminate  which  is 
mid-plane  symmetric  initially  may  behave  in  an  unsymmetrical  fashion 
when  subjected  to  a nonuniform  moisture/ temperature  environment. 

This  is  because  of  different  moisture  contents  and/or  temperature 
through  the  thickness  and  as  a result,  different  material  properties 
for  each  layer.  Thus,  not  only  the  flexural  properties  are  changed 
because  of  the  inclusion  of  moisture/temperature  effects,  but  the 
laminate  also  tends  to  behave  in  an  unsymmetric  fashion  and  thus 
increasing  the  coupling  effect.  Of  course,  since  the  properties  of 
the  lamina  are  functions  of  the  moisture  content  and  the  moisture 
content  is  a function  of  the  thickness  coordinate,  appropriate 
shift  factors  will  be  utilized  from  existing  experimental  data  to 
determine  the  laminate  properties  as  a functipn  of  temperature  and 
moisture  distribution. 

A computer  program  for  numerical  solution  of  the  moisture 
diffusion  equation  with  temperature/moisture  dependent  dif fusivities 
which  vary  across  the  thickness  as  well  as  change  with  time  has  been 
developed.  A finite  difference  scheme  in  space  coordinate  and  a 
predictor  corrector  technique  in  time  variable  are  employed  for 
this  purpose.  The  approach  is  outlined  in  Appendix  D.  Sharp  dis- 
continuity in  dif fusivities  due  to  the  presence  of  layers  of  dif- 
ferent materials  is  also  considered.  For  this  purpose,  moisture 
content  in  the  two  materials  at  the  common  boundary,  as  well  as 
moisture  transfer  across  the  interface,  are  assumed  to  be  equal. 

The  approach  discussed  here  (details  of  which  are  given  in  Appen- 
dix D)  has  been  used  to  solve  the  moisture  diffusion  problem  in  a 
coated  composite.  The  results  are  compared  in  figure  26  with 
those  obtained  in  reference  32.  It  should  be  noted  that  the 


difference  between  the  two  solutions  is  not  significant  and  appears 
to  be  entirely  due  to  the  fact  that  in  reference  32  the  moisture 
concentrations  in  the  two  materials  at  an  interface  are  not  assumed 
to  be  the  same  whereas  they  are  considered  the  same  in  the  present 
study.  The  procedure  outlined  here  can  be  easily  modified  to 
account  for  different  concentrations  in  the  two  materials  as  con- 
sidered in  reference  32.  No  attempt  is  made  in  this  study  to 
obtain  solution  of  the  nonlinear  diffusion  problem  in  the  case  of 
moisture-dependent  dif fusivities  because  data  describing  depen- 
dence of  diffusivity  on  moisture  concentration  are  scarce.  How- 
ever, the  approach  presented  is  a general  one  and  can  yield  solu- 
tion of  the  nonlinear  problem  without  any  more  effort  than  what  is 
needed  in  the  linear  case  as  appropriate  input  data  become  avail- 
able. 
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PANEL  FLUTTER  OF  VISCOELASTIC  COMPOSITE  LAMINATES 


The  capability  to  assess  the  dynamic  response  of  composite 
materials  becomes  increasingly  important  as  the  number  of  their 
applications  in  aerospace  structures  increases.  A number  of  inves- 
tigations, as  noted  earlier,  have  been  performed  in  this  regard  to 
assess  the  influence  of  the  geometric  and  material  variables  upon 
the  vibrational  frequencies  and  dispersion  characteristics.  These 
studies  are  an  essential  first  step  in  assessing  the  response  of 
a structure  subjected  to  a more  complex  transient  and  aerodynamic 
load . 

A direct  and  important  consequence  of  aerodynamic  loading  is 
panel  flutter.  The  problem  of  panel  flutter  (as  distinct  from 
the  more  conventional  lifting  surface  flutter)  can  be  stated  quite 
simply.  A fluid  stream  passes  over  the  surface  of  a thin  plate 
(br  shell)  as  illustrated  in  figure  27.  The  objective  is  to  deter- 
mine the  stability  or  instability  of  the  fluid  structural  system. 

In  particular,  one  is  see)cing  to  establish  the  stability  or  flutter 
boundaries  separating  stable  and  unstable  motion.  If  the  system 
is  unstable,  the  nature  of  the  unstable  solution  can  be  investigated 
further  by  considering  nonlinear  oscillations.  On  the  other  hand, 
if  the  system  is  stable,  the  structure  merely  responds  to  fluctua- 
tions in  the  fluid  stream.  A schematic  of  structural  response  is 
given  in  figure  28 (from  ref.  38)  where  the  plate  deflection  to 
the  plate  thic)cness  ratio  is  plotted  versus  the  dynamic  pressure 
of  the  flow.  Belov;  a certain  value  of  the  dynamic  pressure,  the 
predominant  phenomenon  is  that  of  noise,  while  above  a certain 
dynamic  pressure,  the  amplitudes  increase  significantly  and  that 
is  termed  flutter.  Also,  panel  flutter  is  a supersonic  phenomenon. 
As  is  obvious  from  figure  28  for  larger  values  of  the  dynamic 
pressure,  geometric  nonlinearities  influence  the  problem.  Physic- 
ally, the  nonlinearity  results  from  the  tension  induced  in  the 
plate  by  bending  deformations.  This  tension  increases  with 
increasing  plate  deflection  and  thereby  limits  the  unstable  motion 
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to  a finite  amplitude.  This  limitation  often  causes  the  modes  of 
structural  failure  to  be  those  peculiar  to  fatigue  rather  than 
explosive  fracture  of  the  skin  surface.  A linear  model  (like  the 
one  considered  here)  vjhich  does  not  include  the  effects  of  in- 
plane deformations  of  the  plate  would  of  course  predict  infinite 
deflections  at  the  critical  dynamic  pressure. 

Considerable  work  has  been  done  in  the  problem  of  panel 
flutter  over  the  past  two  decades  (see,  for  example,  refs.  34- 
38).  Most  analyses  can  be  placed  in  one  of  four  categories  based 
on  the  structural  and  aerodynamic  theories  considered;  (i)  linear 
structural  theory;  quasi-steady  aerodynamic  theory;  (ii)  linear 
structural  theory;  inviscid,  potential  aerodynamic  theory; 

(iii)  nonlinear  structural  theory;  quasi-steady  aerodynamic  theory; 
and  finally  (iv)  nonlinear  structural  theory;  inviscid,  potential 
aerodynamic  theory.  As  one  would  expect,  the  significant  body  of 
the  literature  is  in  the  first  category  because  it  is  the  simplest. 
Unfortunately,  this  analysis  has  two  major  weaknesses:  (i)  it 
does  not  account  for  structural  nonlinearities;  hence  it  can  only 
determine  the  flutter  boundary  and  can  give  no  information  about 
the  flutter  oscillation  itself;  and  (ii)  the  use  of  quasi-steady 
aerodynamics  neglects  the  three-dimensionality  and  unsteadiness 
or  memory  of  the  flow.  Hence  it  cannot  be  used  for  Mach  number 
close  to  1.  However,  it  is  known  that  panel  flutter  occurs  in 
the  vicinity  of  M ^1.  Since  it  is  not  the  intent  of  the  present 
effort  to  investigate  the  effect  of  geometric  nonlinearities  as 
well  as  the  appropriateness  of  the  various  aerodynamic  theories, 
only  the  linear  structural  theory  and  the  quasi-steady  aerodynamic 
theory  of  the  first  category  would  be  utilized  for  the  analysis. 
Specifically,  the  linear  piston  theory  aerodynamics  (ref.  39) 
which  includes  aerodynamic  damping  will  be  used  in  the  proposed 
investigation.  This  theory  proposes  that  the  aerodynamic  loading 
on  the  panel  is  proportional  to  the  instantaneous  slope  of  the 
panel  (which  also  suggests  that  panel  flutter  is  a nonconservative 
stability  problem) . 


The  critical  dynamic  pressure  will  be  a function  of  the 
composite  laminate  panel  constitutive  relations  (which  consider 
dynamic  viscoelastic  behavior  of  the  laminate  as  well  as  the 
bending-extension  and  bending-twisting  coupling) . Other  factors 
which  affect  the  critical  pressure  are  the  flow  angle,  and  also 
the  presence  of  in-plane  compressive  forces  (prestress) . In 
the  presence  of  coupling  in  laminates,  the  in-plane  forces  will 
assume  particular  significance  because  their  destabilizing  effect 
will  increase.  The  panel  geometry  (plate  aspect  ratio)  and  boundary 
conditions  will  also  have  noticeable  influence  on  the  flutter 
parameter.  However,  in  the  present  study,  the  primary  objective 
is  to  assess  the  influence  of  some  of  the  unique  characteristics 
of  composite  laminates  on  panel  flutter.  These  characteristics 
are  damping,  coupling,  and  shear  deformation  and  environmental 
effects.  Effects  of  these  characteristics  will  be  investigated 
in  detail  and  the  importance  of  other  analysis  variables  will  be 
studied  in  a cursory  fashion. 

Damping  effects  can  be  considered  for  the  problem  of  panel 
flutter  by  the  inclusion  of  aerodynamic  and  structural  damping 
(for  example,  refs.  40  and  41).  In  the  former  case,  damping  is 
proportional  to  velocity,  while  in  the  latter  case,  damping  is 
incorporated  through  the  complex  moduli.  Inclusion  of  damping 
is  important  for  the  present  problem  because  of  the  following 
reasons : 

(i)  Damping  is  an  inherent  property  of  composite  laminates 
and  as  such  cannot  be  ignored  for  structural  dynamics 
problems . 

(ii)  Though  the  velocity  dependent  damping  increases  the 
flutter  dynamic  pressure,  it  has  been  observed 
that  damping  in  the  form  of  viscoelastic  complex 
moduli  can  frequently  have  a destabilizing  effect. 

(iii)  A theory  which  neglects  damping  predicts  zero  dynanic 
pressure  for  flutter  whenever  the  mid-plane  load  has 
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caused  two  panel  vibration  frequencies  to  coalesce 
(zero  dynamic  pressure  flutter  points) . 

If  one  considers  aerodynamic  and  structural  damping 
where  the  effect  of  the  latter  is  included  by  modifying  both  the 
membrane  and  bending  terms  with  a linear  damping  coefficient,  the 
same  anomolous  behavior  as  in  (iii)  is  predicted.  In  reference  41 
it  was  suggested  that  structural  damping  should  modify  only  those 
terms  in  the  equations  associated  with  bending.  Such  a procedure 
removes  the  physically  untenable  effects  which  render  the  flutter 
parameter  insensitive  at  the  transition  point  to  the  panel  aspect 
ratio  and  boundary  condition.  With  reference  to  the  complex 
properties  of  composite  laminates,  such  a procedure  for  consider- 
ing only  flexural  damping  only  cannot  be  always  true  because  of 
the  possibility  of  extension-bending  coupling. 

The  flutter  of  sandwich  and  orthotropic  panels  has  also  been 
the  subject  of  numerous  investigations  (for  example,  refs.  42-44). 
However,  it  is  only  recently  that  some  effort  has  been  directed 
towards  the  flutter  analysis  of  composite  laminate  plates.  The 
behavior  of  composite  laminates  is  different  from  orthotropic 
panels  (which,  incidentally,  are  a special  case  of  a general 
composite  laminate)  because  of  the  coupling  effects.  Results  in 
references  45  and  46  indicate  that  the  ply  orientations  (and, 
hence,  coupling)  and  the  anisotropy  ratios  strongly  influence 
the  dynamic  instability. 

The  effect  of  shear  flexibility  and  rotary  inertia  has  been 
considered  primarily  for  sandwich  panels  because  of  the  finite 
core  shear  stiffness  (see,  for  example,  ref.  43).  Inclusion  of 
both  shear  flexibility  and  rotary  inertia,  while  still  affecting 
the  flexural  motion,  admits  two  additional  degrees  of  freedom. 

They  are  termed  as  the  thic)cness  shear  and  thic)cness  twist  modes. 
In  general,  both  these  motions  produce  transverse  deflections 
and  thus  influence  the  aerodynamic  loading.  It  was  observed  in 
reference  4 3 that  these  two  modes  mar)cedly  change  the  frequency 


-20- 


coalescence  behavior  and  failure  to  account  for  them  can  lead  to 
significant  overestimates  of  flutter  dynamic  pressure  values. 
Once  again,  for  composite  laminates,  the  "soft"  shear  moduli 
increase  the  importance  of  shear  deformation.  Additionally, 
both  the  shear  deformation  and  material  damping  are  matrix 
dominated.  Hence  their  combined  effect  can  affect  panel  flutter 
significantly . 

In  the  present  investigation,  a "refined"  shear  deformation 
theory  is  utilized  to  characterize  the  composite  laminate  panel 
response.  The  shear  correction  factors  in  this  theory  are 
obtained  for  the  specific  laminate  under  consideration. 

PANEL  FLUTTER  OF  SYMMETRIC  LAMINATES  IN  CYLINDRICAL  BENDING 

The  emphasis  in  the  present  study  is  on  the  composite 
material  parameters  rather  than  on  the  geometric  and  aerodynamic 
aspects  of  the  panel  flutter  problem.  Hence  the  logical  choice 
for  the  panel  geometry  is  the  one  corresponding  to  a cylindrical 
bending  deformation  (infinitely  long  in  the  y direction)  with 
simply  supported  edges.  The  procedure  for  flutter  analysis  of 
such  panels  is  outlined  in  Appendix  C. 

Given  the  frequency-dependent  complex  moduli  of  the  con- 
stituents (fiber  and  matrix)  at  different  temperatures  and  for 
various  moisture  concentrations,  the  lamina  properties  can  be 
computed  by  the  use  of  methods  outlined  in  Appendix  A giving  due 
consideration  to  the  actual  temperature  and  moisture  distribu- 
tion in  the  laminate.  Complex  moduli  of  each  lamina  for  various 
frequencies  are  calculated  for  the  actual  temperature  and  mois- 
ture content  at  the  lamina  location  by  the  use  of  linear  inter- 
polation. One  can  then  obtain  the  laminate  properties  and  shear 
correction  factors  based  on  the  procedure  outlined  in  Appendices 
A and  B.  Formulation  of  the  flutter  determinant  for  a general 
laminated  plate  by  the  use  of  the  Galerltin  procedure  is  dis- 
cussed in  Appendix  C.  In  general,  it  is  necessary  to  evaluate 
eigenvalues  of  a matrix  of  size  5N  x 5N,  where  N is  the  number 


of  terms  retained  in  each  of  the  infinite  series  expressions  chosen 
for  the  five  displacement  variables.  A technique  for  reducing 
the  size  of  the  matrix  by  omitting  some  of  the  modes  which  are 
unimportant  is  also  discussed  in  the  Appendix.  Since  the  lami- 
nate properties  are  dependent  on  frequency,  the  eigenvalues  must 
be  determined  by  the  use  of  trial  and  error  procedure.  At  first 
the  eigenvalues  for  X = 0 corresponding  to  each  N are  determined. 
Attention  is  then  restricted  to  frequencies  less  than  or  equal 
to  the  highest  flexural  frequency.  This  frequency  range  is  then 
divided  into  several  parts.  Average  properties  in  each  part  are 
then  used  to  determine  the  points  at  which  imaginary  parts  of 
eigenvalues  change  sign  as  X is  gradually  increased.  Properties 
at  these  points  (corresponding  to  the  real  part  of  co  ) are  then 
used  to  recalculate  the  eigenvalue  and  the  flutter  parameter. 

This  procedure  is  repeated  to  achieve  the  desired  accuracy. 
Convergence  of  any  trial  and  error  procedure  is  strongly  dependent 
on  the  starting  point  as  well  as  the  increments  by  which  the 
parameter  to  be  determined  is  changed.  The  procedure  used  here 
has  the  same  limitation  and,  therefore,  in  certain  cases,  it  had 
to  be  repeated  to  obtain  the  results  sought. 

An  alternate  approach  to  determine  the  flutter  parameter  for 

(0/90)  laminates  is  also  presented  in  Appendix  C.  In  this 

approach,  the  problem  is  reduced  to  the  calculation  of  the 

determinant  of  a matrix  of  size  N x N for  a range  of  values  of 

X and  real  part  of  co.  Denoting  the  real  and  imaginary  parts  of 

these  determinants  as  A and  A , the  value  of  X can  be  obtained 

R I cr 

from  the  condition  A^^  = A^  = 0.  In  the  initial  stages  of  this 
work,  this  approach  was  used  to  obtain  the  required  value  of  X 

cr 

and  the  results  obtained  are  presented  in  the  next  subsection. 

Results  for  (0/90)  Laminates 
Obtained  by  the  Use  of  the  Criterion  A^^  = A^.  = 0 


Results  for  the  following  three  symmetric  laminates,  each 
0.06"  thick,  are  presented  here  (the  fiber  direction  in  the  0" 
layer  is  assumed  to  coincide  with  the  x axis) ; 
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t ■ (i)  Laminate  1,  Boron/epoxy,  [0/90]^,  span  a = 0.6", 

a/H  = 10 

(ii)  Laminate  2,  Boron/epoxy,  [0/90] span  a = 0.6", 
a/H  =10 

(iii)  Laminate  3,  Graphite  (T300/epoxy)  [0/90] 2^, 

' span  a = 3.6",  a/H  = 60 

All  the  layers  in  the  laminates  are  of  equal  thickness.  Proper- 
ties of  the  fibers  and  the  matrix  are  listed  in  table  1.  Both 
boron  and  graphite  fibers  are  assumed  to  be  elastic  with  the 
same  properties  for  all  frequencies.  The  epoxy  matrix  is  elastic 
under  hydrostatic  stress  with  bulk  modulus  K = 0.3824  x 10^  psi 
for  all  frequencies.  Variation  of  Young's  modulus  E with  frequency 
in  the  range  5 < Re  (to)  < 3 x 10  is  taken  from  reference  2. 

Because  of  the  low  value  of  a/H  ratio  for  the  first  two  laminates, 
properties  at  higher  frequencies  are  required  in  the  flutter 
analysis.  They  are  obtained  (see  table  1)  by  linear  extrapolation 
from  the  data  in  reference  2.  The  value  of  the  shear  correction 
factor  kg^  is  chosen  to  be  0.8.  Since  the  shear  moduli  S3  and 
each  of  the  layers  are  of  the  same  order,  this  choice  is 
adequate.  For  each  of  the  laminates,  several  levels  of  damping 
have  been  considered.  The  imaginary  parts  of  E and  v of  the 
epoxy  matrix  in  table  1 for  all  frequencies  are  multiplied  by  a 
factor  depending  on  the  level  of  damping  (as  listed  in  table  2) 
and  the  real  parts  are  left  unaltered. 

The  A-Re(w)  curves  on  which  the  real  and  imaginary  parts  of 

the  determinants  are  zero  are  plotted  in  figures  29  through  37  • 

t -2 

' for  the  three  laminates  under  consideration.  Note  that  m = 

~ 2 - 

P u)  , where  P is  defined  in  equation  (C3)  of  Appendix  C.  Calcula- 
tions for  Graphite/epoxy  laminates  with  constructions  similar 
to  laminates  1 and  2 and  Boron/epoxy  laminates  with  the  construc- 
tion of  laminate  3 have  also  been  made.  Because  of  similar  trends 
in  the  results,  a discussion  pertaining  to  these  additional 
laminates  has  been  omitted  here.  For  the  same  reason,  results 
for  Graphite/epoxy  and  Boron/epoxy  laminates  of  [0/90]^  construc- 
tion and  length  3.6"  (same  as  that  of  laminate  3)  are  not  pre- 
I sented. 
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For  the  first  two  laminates  (a/H  = 10) , shear  deformations 

have  significant  effect  on  the  modes  of  flutter  as  seen  in 

figures  29-35.  On  the  other  hand,  for  laminate  3 (a/H  = 60) 

and  flutter  occurs  by  coalescence  of  the  first  two  modes  which 

are  predominantly  flexural  in  nature  (figs.  36  and37  ).  In  the 

first  two  laminates,  the  A-Re(cij)  plots  (figs.  30  -32,  34  and35  ) 

show  three  points  of  kinetic  instability.  Two  of  these  points 

are  clearly  due  to  coalescence  of  two  consecutive  modes  of  flexural 

vibration.  Table  3 lists  the  parameters  A . (coalescence  of 

_ cr  1 

fourth  and  fifth  mode)  and  A _ (coalescence  of  second  and  third 

cr^: 

mode)  for  kinetic  instability  for  various  levels  of  damping. 

The  first  flexural  mode  does  not  appear  to  coalesce  with  other 
modes  in  the  range  of  A and  Re((ij)  considered,  but  it  appears  that 
there  is  a point  (^^^.3)  which  unstable  motion  is  possible 
when  damping  is  present.  However,  as  A is  increased  further, 
one  finds  another  point  of  cross-over;  i.e.  the  motion  tends  to 
become  stable  again  (see  fig.  32  ) . This  point  of  cross-over 
occurs  at  a value  of  A which  is  high  as  compared  to  A^^^  in  the 
figure.  The  value  of  A^^^  is  very  low  and  this  instability  does 
not  occur  when  there  is  no  damping. 

Similar  results,  obtained  by  the  use  of  the  method,  outlined 
in  Appendix  C,  for  calculating  real  and  imaginary  parts  of  oj  are 
presented  later.  These  results  indicate  that  as  A is  increased 
beyond  the  imaginary  part  of  u)  changes  signs,  indicating 

that  the  motion  becomes  unstable.  However,  the  oJj.  does  not  in- 
crease rapidly  as  A is  increased,  which  means  that  the  instability 
is  not  of  a very  severe  nature. 

For  laminate  3,  the  frequency  at  which  flutter  appears  to 
occur  is  between  the  first  two  frequencies  of  flexural  motion  and 
the  effects  of  sheer  deformation  do  not  seem  to  have  any  significant 
influence  on  flutter  characteristics  (see  figures  36,  37).  When 
damping  is  present,  the  critical  flutter  parameter  is  significantly 
lowered  (see  table  3) . However,  as  A is  increased  beyond  A^^,  the 
curves  A = f (a)„)  for  A = 0 and  A_  = 0 approach  each  other  again, 

K R X 
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indicating  that  the  instability  becomes  less  severe  as  X is  in- 
creased. Beyond  X = 215  the  two  curves  separate  out  again,  in- 
dicating an  onset  of  severe  flutter.  These  characteristics  are 
clearly  demonstrated  in  the  results  presented  in  the  next  sec- 
tion. It  should  be  pointed  out  here  that  possible  existence  of 
a weak  first  mode  flutter  in  low  supersonic,  transonic  range  in 
other  theoretical  predictions  and  experimental  results  similar 
to  those  observed  here  is  discussed  in  reference  38. 

The  results  given  in  this  section  show  that  introduction  of 
little  damping  lowers  the  flutter  load  considerably  from  the  case 
when  there  is  no  damping.  Further  increase  in  damping  causes 
the  flutter  load  to  increase  again  as  observed  in  laminates  1 and 
2.  However,  the  increase  is  not  very  significant.  These  obser- 
vations are  in  agreement  with  well-known  results  of  destabilizing 
effect  of  damping  in  nonconservative  elastic  systems  (see,  for 
example,  refs.  40,  47).  The  destabilizing  effects  observed  here 
as  well  as  the  severity  of  instability  corresponding  to  X^^2 

found  in  laminates  1 and  2 and  X in  laminate  3 have  been  inves- 

cr 

tigated  further  by  employing  techniques  described  in  Appendix  C 
for  determining  real  and  imaginary  parts  of  w.  These  results 
are  presented  in  the  next  sub-section. 

Results  Obtained  from  Calculation  of  Complex  Eigenvalues 

To  study  environmental,  damping,  coupling  and  shear  deforma- 
tion effects  on  strong  as  well  as  weak  instabilities,  complex 
eigenvalues  o)  have  been  calculated  for  the  following  T300/epoxy 
symmetric  laminates.  The  0°  direction  is  assumed  to  coincide  with 
the  x-axis  and  thickness  of  the  laminates  is  equal  to  .06". 

(i)  Laminate  4,  (0/90)2s  laminate  is  reconsidered, 

a/h  = 60,  a = 3.6" 

(ii)  Laminate  5,  (0/90)2^,  a/h  a 10,  a = 0.6" 

(iii)  Larfiinate  6,  (±  45/02/90)  s , a/h  = 40,  a = 2.4" 

(iv)  Laminate  7,  (± 45/O2/9O) s , a/h  = 10,  a = 0.6" 

T300  fiber  properties  used  for  these  calculations  are  the 
same  as  those  listed  in  table  1.  Input  epoxy  matrix  properties 


are  plotted  against  frequency  in  figure  38.  To  calculate  matrix 
properties  at  different  temperature  and  moisture  concentrations, 
use  is  made  of  the  well  known  shift  hypothesis  (references  48, 

49) . Assumed  shift  factors  are  also  shown  in  the  same  figure. 

For  frequencies  in  the  range  10  - 10^  Hertz  in  situ  storage  modulus 
and  loss  tangents  are  assumed  to  be  the  same  as  those  given  in 
reference  2.  For  low  frequencies,  these  values  are  assumed  to 
behave  in  a manner  similar  to  those  of  an  epoxy  resin,  as  re- 
ported in  reference  48  (with  due  consideration  to  shifts) . The 
bulk  modulus  of  the  matrix  is  assumed  to  be  0.488  x 10^  psi 
( from  ref.  2)  which  is  real  and  constant  for  all  temperature 
and  moisture  content.  It  should  be  pointed  out  that  this 
value  is  different  from  the  one  listed  in  table  1,  used  to 
obtain  the  results  in  one  previous  subsection,  and  for  this  reason 
quantitative  results  presented  here  may  differ  a little  from  those 
discussed  earlier.  Fiber  properties  are  assumed  to  be  real  and 
independent  of  temperature  and  moisture  content.  Experimental 
data  on  constituent  properties  at  various  frequencies,  temperature 
and  moisture  content  are  scarce,  and  therefore,  the  approach 
which  is  discussed  above  and  appears  to  be  a consistent  one,  has 
been  used  to  estimate  the  input  data.  The  results  obtained  are 
discussed  and  conclusions  are  drawn  with  due  consideration  to 
the  calculated  laminate  properties  and,  therefore,  effects  of 
other  types  of  input  properties  can  be  estimated  if  the  corres- 
ponding laminate  properties  are  known. 

Real  and  imaginary  parts  of  io  = tu.  -j  Pa^/D^  for  various  values 
of  \ = 2a^q/(S  D^)  are  plotted  in  figures  39-42  for  the  laminates 
4-7,  respectively.  The  values  of  critical  flutter  parameter  are 
also  listed  in  these  figures.  The  solid  lines  represent  the 
solutions  for  the  reference  state,  i.e.  the  laminate  is  at  room 
temperature  (70°F.)  and  is  dry.  The  dotted  lines  are  for  a pre- 
scribed moisture  and  temperature  gradient  as  shown  in  figure  39. 

The  temperatures  on  top  and  bottom  surfaces  are  250®  and  70 ®F 
respectively.  Moisture  concentration  in  the  matrix  materials 
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on  these  surfaces  are  2.57%  and  zero.  Five  terms  of  the  series 
representation  (N  = 5)  were  retained  in  all  the  calculations.  A 
total  of  twenty  modes  to  include  the  effect  of  four  displacement 
variables,  W,  ip  + W,  x,  ip  and  U , were  considered  in  laminate  7, 

X y 0 

where  bending-extension  coupling  effects  due  to  unsymmetric  mois- 
ture, temperature  distribution  should  have  the  largest  influence. 
No  destabilizing  effect  was  noticed.  Extensional  stiffnesses 
appear  to  be  extremely  high  to  have  any  influence  on  flutter 
characteristics.  Therefore,  fifteen  modes  were  retained  for 
laminate  6 and  ten  were  retained  for  laminates  4 and  5,  where 
there  is  no  coupling  between  the  shear  deformations  + W,x 
and  ip  . 

y 

Observations  and  comments  made  in  the  previous  subsection 
also  hold  for  the  results  presented  in  figures  39-42.  Figure  40 
clearly  shows  that  the  imaginary  part  of  to  increases  slowly  and 
then  decreases  again  as  A is  increased  beyond  Instabili- 
ties at  A = and  much  more  severe.  Shear  deforma- 

tion has  a tremendous  influence  on  the  critical  flutter  parameter. 
By  arbitrarily  increasing  the  shear  stiffness  (increasing  the 
shear  correction  factor  by  1000  times) , the  critical  flutter 
parameter  for  laminate  5 is  increased  from  48  to  117  (see  fig. 

40)  and  severe  flutter  then  occurs  by  coalescence  of  modes  1 
and  2 . 

For  laminate  4,  where  shear  deformation  effects  are  neg- 
ligible, A^^  is  significantly  lowered  by  introduction  of  damping. 
Also,  the  imaginary  parts  of  ui  change  sharply  from  positive  to 
negative  at  but  severity  of  flutter  does  not  increase  until 

X is  increased  by  quite  a good  amount  beyond  This  charac- 

teristic has  been  pointed  out  in  the  previous  subsection. 

In  laminate  7,  the  wea)c  instability  sets  in  at  very  low 
value  of  A as  in  the  case  of  laminate  5.  However,  as  A is 
increased  significantly,  ui^  corresponding  to  mode  1 changes  sign 
again  from  negative  to  positive,  but  corresponding  to  mode  2 
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becomes  negative  and  yields  a severe  instability.  On  the  other 
hand,  it  appears  that  inclusion  of  moisture-temperature  effects 
causes  corresponding  to  mode  1 to  increase  in  magnitude  causing 
a strong  flutter  instability.  Near  these  points  of  severe 
instability,  strong  coupling  exists  between  the  mode  shapes  and 
it  is  difficult  to  accurately  determine  which  frequency  corre- 
sponds to  each  mode.  Calculations  were  made  to  keep  track  of 
the  modes  and  associated  frequencies  and  the  results  appear  to 
confirm  the  characteristic  described  above.  Although  this  behavior 
does  not  influence  critical  values  of  flutter  parameters  for  weak 
or  strong  instabilities,  it  is  an  interesting  phenomenon  and 
should  be  investigated  in  detail  in  future  studies. 

Moisture-temperature  gradient  does  not  appear  to  have  much 

influence  on  the  critical  value  of  the  flutter  parameter  for 

laminate  4,  but  increases  the  critical  values  for  laminates  5-7. 

This  beneficial  effect,  however,  is  not  due  to  increased  damping. 

Moisture-temperature  gradient  actually  reduces  the  loss  tangents 

of  bending  as  well  as  shear  stiffness  of  all  laminates  as  is  clear 

from  the  values  of  at  X = 0.  This  is  a consequence  ot  .:he  type 

of  master  curve  assumed  (fig.  38)  and  the  shift  hypothesis.  Loss 

tangents  for  the  epoxy  matrix  reduce  with  decreasing  frequency 

up  to  (1)  » 10.  Increased  temperature  and  moisture  effects  cause 

a shift  towards  the  left  in  figure  38  by  the  shift  hypothesis, 

which  causes  a reduction  in  loss  Langent.  For  the  same  reason, 

the  storage  modulus  increases  due  to  moisture-temperature  effects. 

This  increase  in  storage  modulus  increases  the  shear  stiffness. 

It  appears  that  this  increase  in  shear  stiffness  is  reflected  in 

the  increased  values  of  X in  laminates  5-7  where  shear  deforma- 

cr 

tions  have  significant  influence.  Since  X^^  for  laminate  4 is  not 
strongly  influenced  by  shear  deformations,  no  increase  is  noticed 
in  this  case. 

The  approach  based  on  calculation  of  complex  eigenvalues 
should  yield  values  of  X^^  which  are  more  accurate  than  those 
obtained  in  the  previous  subsection  by  the  use  of  the  condition 
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= 0.  Therefore,  the  effect  of  various  damping  levels 
(table  2)  was  reinvestigated  by  the  use  of  these  refined  calcula- 
tions. Results  are  tabulated  in  table  4.  For  laminate  5, 
and  ^^^2  to  increase  as  damping  level  is  increased  from  ' 

1 to  3.  For  no  damping  there  does  not  appear  to  be  any  point 
of  incipient  flutter  corresponding  to  X^^2  for  cases  with  damping. 
This  point  ^^^.3  is,  however,  not  a point  of  severe  instability. 

For  laminate  4 X^^  sharply  drops  from  the  value  of  340  for  zero 
damping  to  157,  150,  150  for  damping  levels  1,  2,  and  3,  respec- 
tively. The  reason  for  a higher  value  of  X^^  for  level  1 as  com- 
pared to  those  for  levels  2 and  3 is  not  clear.  This  trend  is 
different  from  that  in  other  results  in  tables  3 and  4,  which 
show  some  increase  or  no  change  in  X^^  as  damping  is  increased 
from  level  1 to  3.  It  should  be  noted,  however,  that  the 
imaginary  part  of  X sharply  changes  from  negative  to  positive 
in  the  case  of  laminate  4 (see  fig.  39)  as  it  passes  through 
X and  the  properties  used  for  calculations  are  obtained  by 
linear  interpolation  from  those  given  in  a tabular  form.  In  view 
of  these  facts,  the  small  difference  discussed  above  does  not 
appear  to  be  significant.  In  general,  the  conclusions  which  can 
be  drawn  from  the  present  results  are  the  same  as  those  discussed 
in  the  previous  subsection.  It  should  be  pointed  out,  however,  that 
the  drop  in  X^^  by  introduction  of  small  damping  is  significantly 
higher  when  shear  deformation  effects  are  small  (a/h  is  large) 
as  compared  to  the  cases  when  such  effects  are  dominant  (a/h  is 
small) . 

A simple  way  to  assess  the  severity  of  flutter  is  to 

examine  how  the  magnitude  of  the  imaginary  part  of  u)  increases 

after  it  becomes  negative.  The  imaginary  parts  of  o)  are  plotted 

in  figures  39-42  for  laminates  4-7  which  clearly  show  that  severe 

flutter  does  not  always  set  in  at  the  values  of  X given  in 

cr 

these  figures.  Another  approach  of  judging  whether  the  instability 
is  severe  or  not  is  to  use  the  concept  of  amplitude  ratio  Ap  or 
logarithmic  increment  (refs.  40  and  47).  The  amplitude  ratio 
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A is  defined  as  the  increase  in  amplitude  in  one  cycle,  which  is 
2Tra)  /o) 

equal  to  e -i-  « g^nd  the  logarithmic  increment  is  Log  A . 

S R 

LoggAR  is  plotted  against  A for  laminate  6 in  figure  43,  which 
clearly  shows  that  beyond  A^^  the  amplitude  does  not  increase 
too  rapidly  with  A.  The  values  of  A near  which  amplitudes  are 
expected  to  increase  rapidly  are  higher  than  A^^  and  should  be 
close  to,  but  less  than,  the  critical  flutter  parameter  for  no 
damping.  These  observations  are  in  agreement  with  those  in 
reference  47.  Figure  43  also  shows  the  effect  of  environment 
(moisture-temperature  gradient) , aerodynamic  damping,  and  pre- 
stress. Moisture-temperature  gradient  appears  to  stabilize  the 
system.  It  has  also  been  pointed  out  earlier  that  structural 
damping  effects  are  less  when  moisture-temperature  gradient  is 
present  in  the  laminate,  which  is  a consequence  of  the  master 
curve  of  epoxy  used  in  calculations.  The  stabilization  effect 
is  caused  by  (i)  reduced  structural  damping  and  (ii)  increased 
shear  stiffness,  but  the  increase  in  A^^  is  possibly  entirely 
due  to  increase  in  shear  stiffness. 

Figure  43  also  illustrates  that  A^^  (A  at  loggAR=0 ) is 
significantly  increased  by  the  introduction  of  aerodynamic  damping 
but  the  amplitude  ratio  versus  A plot  does  not  differ  significantly 
from  the  one  with  no  aerodynamic  damping  for  high  values  of  A . 

The  stabilization  effects  of  aerodynamic  damping  have  been 
observed  in  other  studies  (ref.  40).  Effects  of  prestress  have 
been  studied  for  two  values  of  the  ratio  of  prestress  to  static 
buckling  load,  namely  0.5  and  0.8.  A^^  is  not  significantly 

altered  due  to  the  presence  of  prestress,  but  severe  flutter  sets 
in  at  loads  lower  than  that  with  no  prestress.  Effects  of  environ- 
ment, prestress,  and  aerodynamic  damping  on  A^^  are  given  in 
table  5. 

The  following  important  observations  can  be  made  from  the 
results  presented  in  this  and  the  previous  subsection; 
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Small  damping  lowers  the  flutter  load  from  the  value 
for  no  damping.  The  reduction  is  very  significant 
when  shear  deformation  effects  are  small  or  span-to- 
depth  ratio  of  the  panel  is  large.  For  small  span- 
to-depth  ratios,  this  reduction  is  not  that  significant. 
In  general,  higher  damping  increases  the  flutter  load. 
Shear  deformation  effects  have  strong  influence  on 
analytical  estimates  of  the  flutter  characteristics  of 
a laminate. 

When  damping  is  present,  instability  is  not  always 
severe  at  points  of  incipient  flutter.  Values  of  A 
where  instability  is  severe  are  usually  close  to  A 

cr 

for  no  damping. 

For  the  type  of  master  curve  and  shift  hypothesis  used 

in  this  study,  moisture-temperature  gradient  in  a 

laminate  may  increase  the  value  of  A^^  as  well  as  A for 

a prescribed  amplitude  ratio  A . However,  the  results 

R 

are  strongly  dependent  on  the  type  of  master  curve  and 

shift  characteristics  of  the  constitutents . 

Aerodynamic  damping  stabilizes  the  system  and  significant 

increase  is  noticed  in  the  critical  vai.ue  of  A as 

cr 

the  damping  is  increased.  Prestress  destabilizes  the 
system  but  the  effect  on  A^^  is  not  pronounced  for 
values  of  prestress  less  than  the  static  buckling  load. 
Bending-extension  coupling  effects  introduced  due  to  the 
presence  of  moisture-temperature  gradient  do  not  have 
any  significant  influence  on  the  flutter  characteristics 
of  symmetric  laminates. 


CONCLUDING  REMARKS 


The  problem  of  panel  flutter  of  epoxy  matrix  composite 
laminates  is  investigated  here.  Emphasis  is  placed  upon  incor- 
porating into  the  analysis  an  appropriate  constitutive  relation 
for  the  composite  laminate.  This  requires  consideration  of 
coupling  effects  between  extension  and  bending  and  bending  and 
twisting.  It  also  requires  consideration  of  shear  deformation 
effects  as  a function  of  laminate  construction  and  stacking 
sequence,  and  frequency-dependent  dynamic  viscoelastic  behavior 
of  composite  laminates.  In  addition,  an  approach  for  the  deter- 
mination of  the  effects  of  moisture  and  temperature  on  the 
laminate  constitutive  equations  is  outlined.  Finally,  the  panel 
flutter  of  symmetric  laminates  in  cylindrical  bending  is  inves- 
tigated in  order  to  assess  the  effects  of  the  composite  material 
parameters  such  as  laminate  construction,  shear  deformation,  and 
damping  upon  the  critical  flutter  parameter. 

Effects  of  aerodynamic  damping  and  prestress  are  also 
studied  in  a cursory  fashion.  Although  the  results  given  here 
are  for  symmetric  laminates  only,  the  methods  utilized  are 
applicable  to  cylindrical  bending  of  laminates  which  are  not 
necessarily  symmetric. 

The  correspondence  principle  for  vibrations  of  viscoelastic 
composites  is  utilized  for  the  determination  of  lamina/laminate 
complex  properties  from  the  constituent  complex  properties.  It 
is  observed  that  if  the  effective  elastic  moduli  are  of  a form 
such  that  they  are  explicitly  obtained  from  the  constituent  proper- 
ties, the  dynamic  correspondence  principle  is  particularly  simple 
to  apply.  Typical  results  for  lamina/laminates  show  that  their 
matrix-dominated  stiffnesses  are  essentially  governed  by  the  epoxy 
matrix  complex  properties. 

An  analytical/numerical  procedure  to  obtain  the  shear  correc- 
tion factors  associated  with  the  transverse  shear  stiffnesses  is 
developed.  The  approach  consists  of  obtaining  the  cut-off  fre- 
quencies of  propagation  of  plane  waves  in  the  thickness-shear 
mode  for  a laminated  plate  from  a "refined"  shear  deformation 
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theory  as  well  as  "exact  elasticity  solution"  and  then  equating 
them  to  obtain  unknown  values  of  the  shear  correction  factors. 

In  addition,  a method  is  also  formulated  for  evaluating  the  trans- 
verse shear  coupling  term.  The  shear  correction  factors  for  a 
limited  number  of  composite  laminates  are  presented.  It  is  seen 
that  as  the  number  of  layers  in  a laminate  increases^  the  shear 
correction  factors  approach  the  classical  value  of  ^ for  a homo- 
geneous plate,  and  the  magnitude  of  the  shear  coupling  term 
decreases.  The  ratio  r (*^23^^12^  also  influences  the  magnitude 
of  the  shear  correction  factors.  For  values  of  r not  close  to 
unity  and  for  a small  number  of  layers  in  a laminate,  the  correc- 
tion factors  may  differ  significantly  from  the  classical  value. 

The  investigation  of  panel  flutter  characteristics  of 
symmetric  laminates  indicates  that  both  shear  deformation  and 
damping  have  significant  effects.  In  general,  introduction  of 
little  structural  damping  lowers  the  flutter  load  considerably 
from  the  case  when  there  is  no  damping.  Further  increase  in 
damping  causes  the  flutter  load  to  increase.  These  observations 
are  in  agreement  with  the  well-known  results  of  destabilizing 
effect  of  damping  in  nonconservative  elastic  systems.  Introduc- 
tion of  aerodynamic  damping,  on  the  other  hand,  stabilizes  the 
system  and  compressive  prestress  has  a destabilizing  effect.  In 
cases  where  effects  of  shear  deformation  are  small,  i.e.  sp'^n- 
to-depth  ratio  is  large,  the  destabilization  phenomenon  caused 
by  small  structural  dam.ping  discussed  above  is  more  pronounced 
than  that  in  laminates  with  lower  span-to-depth  ratio.  Severe 
instability  does  not  always  set  in  at  points  of  incipient  flutter 
when  damping  is  present.  Effects  of  environment  on  flutter  charac- 
teristics are  strongly  dependent  on  the  type  of  master  curve  and 
shift  hypothesis  used  to  estimate  constituent  properties  for 
various  temperature  and  moisture  concentration  levels.  Experi- 
mental data  on  constituent  properties  are  very  limited  and  only 
a small  number  of  problems  have  been  studied  here.  However,  the 
usefulness  of  the  analytical  techniques  presented  to  study  various 
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environmental  factors  has  been  clearly  demonstrated.  Bending 
extension  coupling  effects  introduced  due  to  nonsymmetric  tempera- 
ture and  moisture  distribution  do  not  appear  to  have  much  influ- 
ence on  flutter  characteristics  of  symmetric  laminates. 

Finally,  a solution  procedure  for  the  moisture  diffusion 
equation  in  the  case  of  temperature/moisture  dependent  variable 
dif fusivities  across  the  thickness  has  been  developed  by  the  use 
of  a finite-difference  scheme  and  a predictor-corrector  technique. 
Such  an  analysis  in  conjunction  with  the  available  data  in  the 
literature  for  moisture/temperature  effects  on  lamina/laminate 
properties  would  provide  the  means  of  assessing  the  effects  of 
environment  on  panel  flutter,  under  transient  as  well  as  steady - 
state  conditions. 


APPENDIX  A. 

EFFECTIVE  COMPLEX  MODULI  OF  LAMINA/LAMINATE 


CORRESPONDENCE  PRINCIPLE  FOR  VIBRATIONS  OF  VISCOELASTIC 
COMPOSITES 

A theory  of  quasi-static  viscoelastic  behavior  of  fiber- 
reinforced  materials  has  been  developed  (see,  for  example,  ref.  1). 
The  basic  tool  in  the  analysis  is  the  correspondence  principle 
which  relates  effective  elastic  moduli  of  composites  to  effective 
relaxation  moduli  and  creep  compliances  of  viscoelastic  composites. 
Subsequently,  it  was  demonstrated  in  reference  1 that  the  effec- 
tive complex  moduli  (characterizing  the  dynamic  behavior  of  visco- 
elastic composites)  can  also  be  related  to  the  effective  elastic 
moduli  of  composites  by  a correspondence  principle.  It  is  worth- 
while here  to  reconstruct  the  logic  developed  in  reference  1 to 
establish  this  dynamic  viscoelastic  correspondence  principle. 

The  effective  relaxation  moduli  of  a statistically  homogeneous 
viscoelastic  composite  may  be  defined  by 
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where  o. . and  e. . are  the  average  stress  and  strain  tensors, 

ID  ID 

respectively,  and  is  the  effective  relaxation  moduli  ten- 

sor. Assume  that  the  average  strain  field  varies  with  time  in 
the  form 


e . . = e . . e 
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where  u)  is  the  frequency  and 
complex  and  a function  of  w 
equation  (Al)  results  in 


is  time  independent  but  may  be 
. Substitution  of  equation  (A2)  in 


-35- 


a . . (t) 
ID 


ia}t 


ID 


T)e  dt, 


If  t - T = t',  one  obtains 
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Equation  (A5)  resembles  the  elastic  stress-strain  law  and  D. .,5 

•k  ^ D ^ 

are  the  effective  complex  moduli.  separated 

into  real  and  imaginary  parts.  Thus 


Dr^,.„  (id))  = of  ..  0 (di)  + iD*  „ (d)) 
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The  real  part  of  equation  (A7)  is  called  the  storage  modulus 
while  the  imaginary  part  is  termed  as  the  loss  modulus.  The 
ratio  of  the  imaginary  part  to  the  real  part  is  defined  as  the 
loss  factor. 

Consider  now  the  Laplace  transform  (LT)  of  the  effective 
relaxation  moduli 


=rjk«<p>  - ’ <''8> 

0 

Then  the  LT  of  equation  (Al)  is  given  by: 
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where 
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^ijk£  transform  domain  (TD)  effective  moduli.  From  equa- 


tions (A5) , (A8) , and  (AlO) , it  can  be  concluded  that 


and 


(All) 


The  static  correspondence  principle  for  viscoelastic  com- 
posites states  that  the  effective  TD  moduli  of  a visco- 
elastic composite  are  obtained  by  the  replacement  of  phase 
elastic  moduli  by  TD  phase  moduli  in  the  expression  for  the  effec- 
tive elastic  moduli  of  a composite  with  identical  phase  geometry. 
In  other  words,  if  the  effective  elastic  moduli  is  defined  as 


®C*  ..  p = ^C*  p {phase  moduli,  phase  geometry) 


(A12) 


then  the  effective  TD  moduli  are 


1,0  (p)  = ..  p{  TD  phase  moduli,  phase  geometry) 


(A13) 


From  equations  (All)  and  (A13) , it  is  evident  that 


D*..p(i(D)  = ®C*  { TD  phase  moduli  (iuj)  , phase  geometry). 
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But  the  TD  phase  moduli  on  the  RHS  of  equation  (A14)  are  merely 
the  same  as  the  phase  complex  moduli  from  equation  (All) . Con- 
sequently, equation  (A14)  is  rewritten  as 


A 0 * i 

Dijk£{i*^)  ~ ^ij)cf  ^ phase  complex  moduli,  phase  geometry; 


(A15) 


Thus,  the  effective  complex  moduli  of  a viscoelastic  composite  are 
obtained  by  replacing  the  phase  elastic  moduli  by  the  phase  com- 
plex moduli  in  the  expressions  for  the  effective  elastic  moduli 
with  identical  phase  geometry.  If  the  effective  elastic  moduli 
are  of  a form  such  that  they  are  explicitly  obtained  from  the 
constituent  properties,  the  above  correspondence  principle  is  par- 
ticularly simple  to  apply.  The  complex  phase  moduli  are  generally 
obtained  as  a function  of  frequency  from  dynamic  tests.  In  what 
follows,  it  will  be  assumed  that  the  complex  phase  moduli  of  fiber 
and  matrix  are  known  for  the  frequency  range  under  consideration 
and  also  for  various  temperature  and  moisture  content. 

LAMINA  COMPLEX  MODULI  FROM  CONSTITUENT  COMPLEX  MODULI 

To  obtain  the  complex  moduli  of  unidirectional  composite 
lamina,  the  correspondence  principle  of  vibration  of  viscoelastic 
composites  as  discussed  above  and  in  reference  1 will  be  utilized. 
If  analytical  expressions  for  effective  elastic  moduli  are  known, 
it  is  possible  to  obtain  the  effective  complex  moduli  in  a rather 
straightforward  manner.  In  general,  closed-form  analytical 
expressions  for  effective  elastic  moduli  are  difficult  to  obtain. 
Several  approximate  expressions  are  available  in  literature  and 
use  of  these  relations  is  likely  to  yield  approximate  expressions 
for  the  complex  moduli.  The  closed-form  expressions  of  the  upper 
bound  results  for  the  effective  elastic  moduli  based  on  the  com- 
posite cylinder  assemblage  model  (refs.  7 and  8)  will  be  used  here. 
It  is  well  known  that  for  a transversely  isotropic  unidirectional 
composite,  the  composite  cylinder  assemblage  model  yields  coin- 
cident upper  and  lower  bounds  for  four  of  the  five  independent 
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moduli.  Also,  the  upper  bound  results  for  the  other  modulus  are 
found  to  be  close  to  experimental  results.  Therefore,  the  results 
given  here  should  yield  good  approximate  results  for  the  complex 
moduli.  It  should  also  be  pointed  out  that  the  expressions  for 
four  of  the  five  independent  moduli  are  correct  within  the  limita- 
tions of  the  composite  cylinder  assemblage  model. 

The  results  obtained  by  the  use  of  the  procedure  described 
above  are  given  below. 
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4K 


K + (1  + 


where  V*"  and  V^  are  the  volume  fractions  of  matrix  and  fibers  in 
the  fiber  bundle,  respectively,  and 
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In  the  above  equations,  E,  v,  G,  and  K are  the  Young's  modulus, 
Poisson's  ratio,  shear  modulus,  and  plane  strain  bulk  modulus, 
respectively.  Superscripts  "f"  and  "m"  refer  to  the  fiber  and 
matrix  material  and  subscripts  "A"  and  "T"  refer  to  the  fiber 
bundle  axial  and  transverse  directions.  Superscript  denotes 


the  effective  fiber  bundle  property.  The  fibers  and  matrix  mate- 
rials are  considered  to  be  transversely  isotropic.  Also,  the 
symbol  over  a quantity  indicates  that  the  quantity  is  complex. 
The  complex  moduli  for  a composite  at  a given  frequency,  tempera- 
ture, and  moisture  content  can  be  obtained  by  using  the  expressions 
given  above. 

LAMINATE  COMPLEX  STIFFNESSES 

Approximate  expressions  for  effective  complex  stiffnesses 
A. . and  D. . for  a laminated  plate  are  obtained  by  the  use  of 
the  displacement  field  in  equation  (B13)  of  Appendix  B,  the 
laminated  plate  theory  of  reference  24,  and  the  correspondence 
principle  outlined  above.  The  expressions  for  these  complex 
stiffnesses  are  given  below: 
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and  (w)  = (i^) 


(i, j = 4,5) 


where 
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Q.  . (uj)  = C.  . (w)  - ;s:: C (to)  (plane  stress  reduced 
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stiffnesses;  i,j  = 1,2, 3,6; 
a = 1,2,6) 


and 


C.j(to)  = frequency-dependent  lamina  complex  stiffness  matrix, 

k. . = shear  correction  factors  for  transverse  shear 

which  are  functions  of  laminate  construction 
and  stacking  sequence 

The  method  of  calculation  of  effective  properties  for  the  laminate 
in  shear  (i,j  = 4,5)  is  discussed  in  Appendix  B. 
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APPENDIX  B. 

DETERMINATION  OF  SHEAR  CORRECTION  FACTORS 
FOR  COMPOSITE  LAMINATES 

One  of  the  objectives  of  the  present  program  is  to  develop 
laminate  constitutive  relations  which  consider  the  effect  of 
laminate  stacking  sequence.  In  a conventional  laminate  analysis, 
the  bending  (D^ ^ , i , j=l , 2 , 6)  and  bending-extension  coupling 
(B^j ,i, j=l,2,6)  stiffnesses  are  functions  of  stacking  sequence 
while  the  extensional  stiffnesses  (A^ ^ , i , j=l , 2 , 6 , 4 and  5)  are 
not.  Hence,  the  stiffnesses  associated  with  the  transverse  shear 
deformation  (A..,  A._  and  A_c)  are  insensitive  to  the  manner  in 
which  the  laminae  are  stacked  together.  In  addition,  the  shear 
correction  factor,  k,  used  with  the  A^^,  A^^  and  A^^  stiffnesses 
is  generally  a constant  (tt^/12  or  5/6)  . An  approach  to  consider 
the  stacking  sequence  effect  is  to  formulate  a higher  order  plate 
theory.  But  this  would  be  inappropriate  and  cumbersome  for  use  in 
structural  dynamics  problems  such  as  panel  flutter.  An  alternate, 
and  perhaps  more  convenient,  approach  is  to  determine  the  shear 
correction  factors  as  a function  of  laminate  construction  and 
stacking  sequence  by  using  a refined  or  "exact  elasticity"  theory 
and  to  utilize  these  shear  correction  factors  in  conjunction  with 
a simple  shear  deformation  laminate  plate  theory  of  the  type  in 
references  12-14  and  23.  The  analytical/numerical  procedure  to 
obtain  the  correction  factors  (associated  with  the  transverse  shear 
stiffnesses  A^^,  A^^  and  A^g)  is  described  subsequently.  The  basic 
approach  is  to  obtain  the  cut-off  (infinite  wave  length)  frequen- 
cies for  propagation  of  plane  waves  in  the  thickness-shear  mode 
for  a laminated  plate,  both  from  the  shear-deformation  theory  and 
the  "exact  elasticity  solution"  and  then  equate  them  to  obtain 
the  unknown  values  of  the  shear  correction  factors  k^^  and  k^g. 

In  addition,  a method  has  been  formulated  for  evaluating  the 
transverse  shear  coupling  term 
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THICKNESS  SHEAR  MODE  CUT-OFF  FREQUENCIES  FROM  THE  ELASTICITY 
SOLUTION 

Consider  a laminated  plate  consisting  of  m layers  of  total 
thickness  H,  as  shown  in  figure  2.  Assume,  . for  example,  the 
case  of  plane  waves  propagating  in  the  x-direction.  Then  the 
displacements  are  of  the  form 

u = u(x,z),  V = v(x,z),  w = w(x,z)  (Bl) 


The  constitutive  relations  for  each  orthotropic  layer  (at  a 
fiber  orientation  6 with  respect  to  the  x-axis)  in  the  laminate 
coordinate  system  xyz  are 
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where  the  strain  components  are  the  engineering  strains.  From 
the  linear  stress-strain  relationships  and  the  equations  of  motion 
from  the  theory  of  elasticity,  the  governing  field  equations  for 
each  layer  are 
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where  a comma  followed  by  a subscript  indicates  a differentiation 
and  p is  the  mass  density.  The  boundary  conditions  on  the  laminate 
surfaces  (z  = ±H/2)  are: 
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The  continuity  of  displacements  and  stresses  at  the  interface 

requires  that  a , i , T , u,  v,  and  w be  continuous  at  the 
Z Z X z ^ ^ 

interfaces . 

A solution  for  the  displacements  in  each  layer  for  the  problem 
of  plane  wave  propagation  can  be  assumed  to  be 
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where  f = tt/Iio  and  A is  the  wave  length.  For  infinitely  large 
wave  lengths  and  for  the  thickness  shear  mode  of  wave  propagation 
(w  = 0),  equations  (B3)  reduce  to 
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For  the  solution  of  equation  (B7),  only  the  boundary  conditions 
in  equation  (B4)  need  to  be  satisfied,  in  addition  to  the  con- 
tinuity of  ^yz'  ^ layer  interfaces.  It  should 

be  noted  that  equation  (B7)  holds  for  any  direction  of  wave 
propagation. 

After  solving  for  U(z)  and  V(z)  in  equations  (B7),  a relation- 
ship can  be  established  between  the  interface  tractions  (t  and 
Ty^)  and  displacements  (U  and  V)  with  reference  to  a local  car- 
tesian coordinate  system  xyz,  ^(see  fig. 2).  z.  „ is  defined  such 
that  the  layer  thickness  h is  bounded  by  the  planes  z^  = ±h/2. 

The  rexationship  between  the  surface  tractions  and  the  displace- 
ments is 


+ 

T 

u 

XZ 

+ 

+ 

T 

V 

yz 

T 

2(l-d^d2) 

u” 

XZ 

T 

v" 

yz 

!_  _ 

where  the  + and  - signs  indicate  the  quantities  at  Zj  = h/2 
and  Z£  = -h/2,  respectively.  The  elements  of  the  matrix  [K] 

are  listed  below: 


k 

k 

k 

k 

k 

k 


11 

12 

13 

14 

22 

24 


■33 


C^5  (Sj_-tj_)  + C^4d2(s2-t2) 


^21  “ ^34  " ^^43  ""  ^55*^1  ^^1”®1^  ^44^2^^2"®2^ 


■31 


-C-5(t^+Si)  - C;4d2(t2^S2) 


‘41 


*^23  “ ^32  “ ^55*^1  ^^l'*'®!^  ^44*^2  ^*^2‘'’®2  ^ 


44 


^55*^1  ^®l”^l^  + 


(B9) 
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■ • 

where  , t2»  and  S2  are  obtained  from 

^^^44^55“^45  ^^1,2  ^ ^ ^^55'''^44  ^ +^^^55"^44^  ■'’  ^*^45  HBIO) 


^1,2 

= ^45^1, 

2/  (P-C55^44Yi,2) 

^55  = 

^1 ^^55 

^ ‘^2^45) 

<4  = 

^ (c 

2 ^‘"44 

^ ^45^ 

(Bll 

t . = tan ( Y . ^h/2 ) 

^ ^ ; i = 1,2  (B12) 

= cot  ( Y^‘*ih/2 ) 

2 

YT  are  the  characteristic  roots  for  the  differential  equations 

1 . . 2 . 

(B7) . The  expressions  in  equations  (B12)  are  valid  if  Y^(i=l»2) 

as  evaluated  from  equation  (BIO)  are  positive.  If  one  of  them 

is  negative,  then  Y.  in  equation  (B12)  is  calculated  as  the  square 

^ 2 

root  of  the  absolute  value  of  Y^  and  the  trigonometric  functions 
are  replaced  by  their  hyperbolic  counterparts. 

Once  the  stiffness  matrix  [K]  for  each  layer  is  obtained, 
the  global  stiffness  matrix  for  the  laminate  made  of  a finite 
number  of  layers  is  constructed  by  considering  the  free  surface 
and  interface  conditions.  Since  the  laminate  surfaces  are  stress 
free,  the  cut-off  frequencies  w are  obtained  by  equating  the 
determinant  of  the  global  stiffness  matrix  to  zero.  Because  of 
the  complicated  manner  in  which  w enters  in  each  element  of 
the  matrix,  a trial  and  error  procedure  is  employed.  The  local 
and  global  stiffness  matrices  are  to  be  obtained  for  each  trial 
value  of  u).  The  size  of  this  global  stiffness  matrix  is  (2m)x(2m) 
where  (m-1)  is  the  number  of  layers.  The  size  of  the  global 
stiffness  matrix  can  be  reduced  for  the  case  of  mid-plane  symmetric 
laminates  or  for  the  case  when  all  the  laminae  material  axes 
coincide  with  the  laminate  x-y  axes. 
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THICKNESS  SHEAR  MODE  CUT-OFF  FREQUENCIES  FROM  LAMINATED  PLATE 
SHEAR  DEFORMATION  THEORY 

The  shear  deformation  theory  utilized  here  is  the  one 
developed  in  reference  12  for  isotropic,  homogeneous  plates  and 
for  laminated  plates  in  references  13  and' 14.  Referring  again 
to  figure  2,  the  displacement  field  is  assumed  to  be  of  the 
form 


u = u (x,y,t)  + z'|'^(x,y,t) 


V 


= v°(z,y,t) 


+ zi;^y(x,y,t) 


w = w®(x,y,t)  (B13) 

where  u^,  v^,  and  are  displacements  at  the  reference  surface 
(center)  of  the  laminate.  It  is  evident  that  the  local  layer 
rotations  are  assumed  equal  and  "thickness  stretch"  is  ignored. 

The  stress  and  moment  resultants  for  the  laminate,  each  per 
unit  length,  are  defined  in  the  usual  way;  that  is. 


H/2 

o 

h 

II 

f 

xy  X y 

-H/2 

H/2 

‘xy'  ■ 

-H/2 

X 

X y xy  xz'  yz 


(B14) 


By  substituting  the  stress-strain  relationships  of  equation 
(Bl)  and  the  displacement  field  of  equation  (B13)  in  conjunction 
with  the  strain-displacement  relations  of  linar  elasticity, 
in  equation  (B14) , the  plate  constitutive  relations  are  obtained 
as: 
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Aij  = / C^jdz  (i, j = 4,5) , 

-H/2 


Qia  = Sc  = 1.2,3,6;a  =1,2,6) 


{B17) 


k,.,  k^_,  k.-  = shear  correction  factors. 
44  55  45 


In  the  contracted  notation  = o^,  = o^,  = a^, 

o_  = T , and  O-  = T . In  general,  the  shear  correction 
5 xz  6 xy  ’ 

factors  are  assumed  to  be  equal;  that  is,  k^^  = k^^  = k^g  = k. 

Referring  to  equations  (Bl3) , it  is  evident  that  the  trans- 
verse shear  strain  is  constant  across  the  laminate  thickness, 
thus  making  the  interlaminar  shear  stresses  discontinuous  at  the 
interfaces.  This  is  a violation  of  the  continuity  of  stress.  In 
reference  23,  the  shear  stress  is  assumed  continuous  and  constant 


across  the  thickness  and  the  shear  strain  is  allowed  to  be  dis- 
continuous. This  is  achieved  in  the  following  manner.  From 

equation  (B2),  the  shear  strains  y and  y in  a layer  are 

yz  xz  ^ 

related  by  the  relationship  (in  contracted  notation) 


c.  (n») 

= S . . o . 
ID  D 


, i»j  = 4,5 


(B18) 


where  are  the  compliances  of  layer  m.  and  are  contin- 

uous across  the  thickness  of  the  plate,  while  the  strains 
are  discontinuous.  Integration  of  equation  (B18)  yields 


r . = L. .Q.  (i, j = 4,5) 

1 ID  D 


(B19) 


where 


H/2 

r.  - H / Yi  dz, 

-H/2 


L.  . = ^ f sf’l’^dz 
ID 

-H/2 


Qj  = Ha  . . 


can  be  interpreted  as  average  laminate  transverse  shear 
strains.  Equation  (B19),  on  inversion,  results  in 


Q.  = A.  . r. 
D ID  1 


' -1 
where  A. . = L. . 

ID  J-D 


(B20) 


I 
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It  is  now  assumed  that  the  displacement  field  in  equation 
(B13)  represents  the  average  plate  displacements.  Thus 


'04" 
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44  44  45  45 
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_ ^45^45  ^55^55 _ 
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t t I 

v;here  k...  and  k.^  are  the  modified  shear  correction  factors 

44  55  45 

and  are  distinct  from  k^^,  k^^,  and  k^^  in  equation  (B16) . In 
general,  one  may  define  the  following  relationship 
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s: 
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^5 
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(B22) 


(B23) 


The  equations  of  motion  in  terms  of  stress  and  moment  resul- 
tants can  now  be  written  as: 


+ N 

= Pu?,.. 

+ Rlii 

x,x  xy,y 

' tt 

^X,tt 

N + N 

= Pv?.  . 

+ Rip 

xy,x  y,y 
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^y,tt 
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(B24) 

'‘x , X y , y 

M + M 

' tt 

- Q = 1 

+ Ii/; 

X , X xy , y 

X 

'tt  ^X,tt 
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(B25) 


M + M - Q = Rv;.  . + I'l'  . . 

xy,x  y,y  '^y  ' tt  y,tt 


H/2 

2 

where  (P,R,I)  = / p(l,z,z  )dz 

-H/2 


For  the  particular  case  of  plane  wave  propagation  in  the 
X direction,  the  above  equations  of  motions  are  rewritten  in 
terms  of  displacement  variables  (from  ref. 14  ) 


0 

■11^'xx 

''l6'''xx  * 

®ll'^x,xx 

+ B,  = Pu?^^  + Rif) 

16  y,xx  tt  x,tt 
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16'^ 'XX 

^le'^x,  XX 

^ ®66'^y,xx  ~ ^'^'tt  ^ ^’^y,tt 

55<'^x,x  ^ 
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Kyic'l^ 
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®16'''xx  * 
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ll’^XjXX 

* “le^y.xx  - ■'55 '♦x  * "-x’ 

= RuT,.,.  + 

4 5^y  tt  ’^Xjtt 


and 


®16^'xx  ^ ^ee'^'xx  ^ ^16'^x,xx  ^ ^66’^y,xx 


-^5^’^x  '^'x)  - ^44^  = '"'^'tt  ^ ^^,tt 


(B26) 


Solutions  to  these  equations  can  be  written  in  the  following 


form; 
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u = 

y^iu)(t-fx) 
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^giu)(t-fx) 

= 
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and 


w = We 


id)  (t-f  x) 


(B27) 


where  f = tt/Au);  o)  is  the  frequency  and  A the  wavelength. 


For  the  thickness-shear  mode  of  wave  propagation  (A  infinitely 
large  and  w = 0) , equations  (B26)  further  reduce  to 


= 0 

55  X 45  y 

-K.c't'  + - K.  = 0 

45  X 44  y 


(b28) 


v/here  I ' = 1 - R'^/P  . 


Solution  of  equations  (B2P)  yields  the  cut-off  frequencies 


1,2 


21 '01 


’'■'1,2  ~ ^^44  *■^55^  ■ " *^44^  ^^^45^ 


1/2 


(B29) 


Note  again  that  to?  ^ are  independent  of  the  direction  of  wave 

^ f • _C  » * 

oropagation  because  of  the  assumption  that  A approaches  infinity. 
Also,  unless  otherwise  stated,  it  will  be  assumed  that  in  mode  1, 
is  the  predominant  displacement  while  in  mode  2,  4^^  assumes 
that  role. 

If  the  cut-off  frequencies  lo,  _ are  known  from  the  "exact 
elasticity"  solution  as  described  earlier,  then  the  quantities 
and  Kgg  are  obtained  from  equation  (B30);  that  is, 

2 2 r22  22  2n 

55,44  = ^'^‘^1  <‘"l  ■ - ^^^45! 


2K 


provided  that  is  known.  An  approximate  procedure  for  estimating 

is  described  subsequently.  Having  obtained  the  values  of 
Kj-g  and  for  a particular  laminate  in  this  fashion,  the  shear- 

deformation  theory  can  be  used  for  calculating  the  response  of  a 
laminated  plate  under  various  loading  conditions. 
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DETERMINATION  OF  THE  TRANSVERSE  SHEAR  COUPLING  STIFFNESS 


It  is  observed  in  equation  (B30)  that  the  shear  correction 

factors  k,.  and  can  be  obtained  when  the  thickness  shear  cut- 

44  55 

off  frequencies  u?  ^ (from  the  exact  elasticity  solution)  and 

It 

the  coupling  term  ^”^45^45  ^45  ^45  known.  By  using 

the  approach  in  reference  16  , the  term  can  be  obtained  to  be 

K45  = <K44l'55)^^\5-  (B31 

1 

Since  the  stiffness  A^^  is  known,  both  and  can  now  be  | 

determined  from  equation  (B30) . Based  on  the  displacement  field  | 

in  equation  (B13) , ^ is  zero  for  orthotropic,  symmetric  and 

[+0]  laminates  (+0  laminates  having  equal  percentages  of  +0  and 

-0  layers) . Thus,  the  coupling  stiffness  also  vanishes  for  i 

these  laminates.  In  reality,  however,  the  transverse  shear  force  j 

will  induce  both  Y (z)  and  v (z)  shear  strains,  although  1 

“x  xz  yz  I 

the  magnitude  of  the  latter  will  be  smaller  than  that  of  the 
former.  The  actual  coupling  is  expected  to  reduce  for  laminates 
which  are  more  homogeneous  = 0 for  homogeneous,  orthotropic 

laminates) . In  the  present  investigation,  the  magnitude  of  this 
coupling  is  ascertained  from  the  procedure  outlined  below. 

2 

The  two  lowest  values  of  the  cut-off  frequencies  _ are 
obtained  initially  from  the  exact  elasticity  solution.  From  the  ' 

eigenfunctions  corresponding  to  each  of  these  two  frequencies, 
the  following  integrals  are  evaluated  for  each  layer:  1 

[C55(U'"  - U”)  + - V")] 


[C^5{U'^  - U")  + - V')]  (B32) 
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h/2 

/ T dz 
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h/2 
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-h/2 


yz 


These  integrals  are  summed  for  all  the  m layers  of  the  laminate 

i.u)^ 

giving  the  shear  forces  Q e and  Q, 


1,2 


1,2 

modes.  The  average  shear  strains  '1'  and  '1' 

^1,2  ^1,2 

estimated  by  utilizing  equation  (B28).  They  are: 


for  the  two 
are  then 


= Q /I'OJ? 
Xi  X.  1 


= Q /I'oj? 

Yi  y/  1 


i = 1,2 


(B33) 


Hiving  obtained  4'  and  T , equations  (B  22i  are  rewritten  as 

^i  ^i 


55  X.  54  y.  ^x. 

1 1 i 


+ K..*'F  = Q 

45  X.  44  y.  y. 
1 1 


i = 1,2 


(b34) 


The  unknowns  are  K^c*,  K..*,  K,c*,  and  K^,*. 

55  44  ' 45  54 

The  coupling  term  to  be  used  in  equation  (b30)  is  taken  to 
be 


(B35) 


Thus,  the  shear  correction  factors  can  be  obtained  either  from 
equations  (B35)  and  (B30)  or  directly  from  equations  (B34) . 

The  differences  between  the  shear  correction  factors  determined 
from  both  the  approaches  are  insignificant.  Hence,  the  procedure 
for  the  determination  of  in  equation  (B35)  appears  to  be 

reasonable. 

It  should  be  noted  that  in  a few  cases  when  K..*  and  K^t:* 

44  55 

are  almost  the  same,  substitution  of  from  equation  (B35) 
in  (b30)  yields  a very  small  negative  number  under  the  square 


-54- 


root  sign  in  the  right-hand  side  of  equation  (B30).  This  dis- 
crepancy has  been  neglected  in  the  results  reported  herein.  Also, 
for  laminates  with  = 0,  use  of  equation  (B23)  yields  = 0. 

Use  of  a zero  value  for  in  equation  (B30)  yields  and 

Kgg,  which  are  significantly  different  from  and  as 

obtained  from  equations  (B34).  This  is  due  to  the  limitation 
of  the  plate  theory  (which  yields  = 0)  and  can  be  avoided  by 
techniques  described  in  this  section. 


APPENDIX  C. 

PANEL  FLUTTER  ANALYSIS  OF  VISCOELASTIC  LAMINATED  PLATES 

IN  CYLINDRICAL  BENDING 

An  important  objective  of  this  study  is  to  investigate  the 
effects  of  environment,  damping  (dynamic  viscoelastic  properties) 
and  coupling  properties,  and  shear  deformation  in  composite 
laminates  on  panel  flutter.  The  panel  flutter  problem  is  a com- 
plex one  and  is  influenced  by  various  factors  such  as  type  of 
aerodynamic  theory  used,  aspect  ratio  of  the  panel  and  boundary 
conditions,  flow  angle,  prestress,  etc.  For  panels  made  of  com- 
posite laminates,  the  variables  in  the  analysis  are  increased 
significantly.  Therefore,  the  emphasis  here  is  on  the  material 
parameters  rather  than  on  the  aerodynamic  or  geometric  aspects 
of  the  problem.  Hence,  it  is  logical  to  consider  a simple  geometry 
and  an  aerodynamic  theory  which  is  amenable  to  a straightforward 
solution  procedure.  To  that  end,  the  problem  of  panel  flutter 
in  cylindrical  bending  is  investigated  here.  The  aerodynamic 
loading  and  damping  is  determined  from  linear  piston  theory. 

A shear  deformation  laminated  plate  theory  which  considers  unique 
shear  correction  factors  for  a given  laminate  construction  and 
stacking  sequence  (see  Appendix  B)  is  employed  and  material  damp- 
ing is  included  by  considering  complex  lamina  stiffnesses  (see 
Appendix  A) . The  effects  of  these  material  parameters  on  the 
panel  flutter  characteristics  (critical  aerodynamic  pressure  or 
speed)  can  then  be  evaluated. 

GOVERNING  EQUATIONS  AND  METHOD  OF  ANALYSIS 

The  panel  flutter  problem  is  illustrated  in  figure  21,  which 
shows  supersonic  flow  over  a flat  plate  at  an  angle  A with  the 


x-axis.  The  linear,  small-deflection  shear-deformation  laminate 
theory  in  conjunction  with  frequency  dependent  complex  laminate 
properties  is  used  to  characterize  the  panel  response.  Since  the 
panel  is  constrained  to  deform  in  cylindrical  bending,  the  analysis 
variables  are  independent  of  the  y coordinate.  Following  refer- 


ence 14 
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governing  equations  for  panel  flutter  can  be  written  as: 

^ 0 ^ 0 ^ ^ 0 

''u“-xx  ^16'''xx  * ®ll*x,xx  * he^y.xx  ' ‘’“'tt  ''♦x.tt 

^16^'xx  ^ee'^'xx  ®16’^x,xx  ®66'^y,xx  “ ^'^'tt  * 

2c 
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cc('P  + Wf  ) + K.c'/'  = COS  A w, 

55  ^x,x  XX  45’^y,x  B 


+ P cw  + N w,  + Pw 
a 't  XX  XX 
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®ll“°xx  ^le'^'xx  °ll’^x,xx  °16'*'y,xx  ~ ^55^'*'x  ”'x^ 

0 

- K45«<^y  = Ru,^^  + 

®16“'xx  ■*■  ^ee'^'xx  ■*■  ^16'^x,xx  °66’^y,xx  " ^45  ^"^x  ”'x^ 


-K  tlJ  =Rv.  +14/ 

44^y  'tt  y,tt 


where 


p = free  stream  air  density, 
c = free  stream  speed  of  sound, 
q = dynamic  pressure  of  the  airstream, 
B = compressibility  factor  = \/m^  - 1 

M = Mach  number 


= panel  prestress 

XX 

P,R,I=  inertial  terms  defined  in  equation  (B25) 

and 

K.  . ,'b.  . /d.  .(i,j  = 1,2,6)  and  “k^.  . (i,  j = 4,5)  are  the 

i j i j i j ^ w w • ^ j w ^ w 

complex,  frequency-dependent  laminate  stiffnesses  derived 
in  Appendices  A and  B. 

The  aerodynamic  loading  is  obtained  from  the  linear  piston 
theory  (ref.  39) . The  piston  theory  gives  a simple  expression 


(Cl) 


) 
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for  the  aerodynamic  loading  and  has  been  shown  to  be  applicable 
for  Mach  nmnbers  greater  than  1.6  (ref.  38).  The  expression 

^ cosAw,^  is  the  aerodynamic  pressure  term,  while  the  term 
P cw,  represents  aerodynamic  damping.  It  is  also  assumed  that 

cl  U 

the  prestress  is  unaffected  by  the  displacements  in  the  plate. 

The  boundary  conditions  at  x = 0,a  (fig.  21)  for  a simple 
panel  are 


w = 0 , 


and 
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In  what  follows,  the  following  nondimensional  quantities 
will  be  utilized: 


X,  Uq,  Vq,W  = (X,  Uq,  Vq,  w)/a 
Aij,  D^.,  K.j  - 


P,  R,  i = (a^P,  a^R,  a^I)/D, 


, 4 

T 2a  q cos  A „ _ ® j « 
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2„0 
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and 


Dq  = a real  constant  with  the  same  unit  as  ^^i*  con- 

venience, a number  close  to  the  real  part  of  the  bending 
stiffness  corresponding  to  the  first  frequency  in  the 
chosen  frequency  range  may  be  selected. 


The  governing  differential  equations  in  terms  of  displacements 
Uq  , ^ rotations  same  as  equation  (Cl) 

with  all  the  quantities  replaced  by  their  nondimensional  counter- 
parts in  equations  (C3) . 

The  Galerkin  method  will  be  used  to  solve  the  differential 
equations  (Cl) . In  this  particular  application  of  the  method,  a 
series  solution  for  the  dependent  variables  is  chosen  which  satisfies 
all  the  boundary  conditions  and  the  coefficients  of  the  series 
are  determined  so  that  each  term  of  the  series  is  orthogonal 
to  the  exact  solution.  The  series  functions  assumed  for  this 
analysis  are 
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00 
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e Z e cos  ex  x 

n n 

n=l 


(C4) 


where  = nir.  The  boundary  conditions  are  identically  satisfied 
by  the  functions  in  equations  (C4)  . Use  of  the  Galerlcin  pro- 
cedure in  conjunction  with  equations  (Cl) , (C2) , and  (C4)  yields 


2 

{ [ ] — ***  j ]}  {A^}  + ~ ^ 1,2,...  “ (C5) 
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If  one  truncates  each  of  the  series  in  equation  (C4)  at 

n N,  equation  (C5)  yields  a set  of  5N  homogeneous  equations. 

For  a nontrivial  solution,  the  determinant  of  the  coefficients 

a , b , c , d , e^  (n  = 1,2...n)  must  equal  zero.  Standard 
n n n n n 

eigenvalue  routines  cannot  be  used  to  calculate  the  eigenvalues 

u)  ^ (for  given  A,  the  flutter  parameter)  since  the  elements  of 
n 

the  matrix  (K^l  are  frequency  dependent.  Hence,  a trial  and 
error  procedure  must  be  employed. 

In  general,  the  frequencies  are  complex;  that  is, 

(i)  = <i)„  + i u)  (C9) 

The  panel  behavior  is  characterized  by  the  variation  of 
(u)_  + i (i)T)with  X.  Since  the  time  dependence  is  defined  by  the 
term  it  is  clear  that  if  the  imaginary  part, 

is  negative,  a divergent  motion  ensues.  If,  on  the  other  hand, 

(Dj  is  positive,  the  result  is  a damped  harmonic  motion.  Thus, 
flutter  is  considered  to  occur  for  the  lowest  value  of  A at 
which  the  imaginary  part  of  the  frequency  changes  sign,  or  in 
other  words,  vanishes.  This  condition  corresponds  to  sustained 
simple  harmonic  motion  and  represents  incipient  flutter. 

For  A = 0,  the  imaginary  parts  of  all  the  frequencies 
are  positive  (damped  motion) . As  A is  increased  to  the  point 
of  flutter,  the  imaginary  part  of  the  one  of  the  frequencies  will 
be  equal  to  zero.  For  a further  increase  in  A,  the  imaginary 
part  will  become  negative.  However,  in  some  cases,  a phenomenon 
termed  here  as  "weak  coalescence"  (of  two  modes)  may  occur. 

In  such  a case,  the  imaginary  part  of  the  frequency  tends  to 
become  zero  again  with  increase  in  A and  careful  consideration 
should  be  given  to  select  the  correct  flutter  parameter  from 
the  point  of  view  of  practical  applications. 

It  has  been  pointed  out  earlier  that  for  numerical  calcula- 
tions, it  is  necessary  to  truncate  the  system  of  equations  (C5) 
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at  n = N and  investigate  the  possibility  of  coalescence  of  any 
two  of  the  5N  modes.  Usually,  the  real  parts  of  frequencies 
corresponding  to  the  modes  dominated  by  the  in-plane  displace- 
ments u® (c  , n = 1,2...N)  and  (e  , n = 1,2...N)  as  well  as  the 
n n 

average  rotations  due  to  shear  + w,  (b  , n = 1,2...N)  and 

ifiy  (d^,  n = 1,2...N)  are  high.  Therefore,  the  possibility  of 

flutter  as  a result  of  coalescence  of  a majority  of  the  modes 

dominated  by  c^,  e , b , and  d (n  = 1,2...N)  with  each  other 
n n n n 

or  with  those  in  which  the  terms  a^  (corresponding  to  the  dis- 
placement w)  are  dominant  can  be  ruled  out  and,  except  for  the 

first  few  of  the  constants  b , c , d , and  e (n  = 1,2,...N), 

n n n n 

all  others  may  be  assumed  to  be  zero  without  any  loss  in 
accuracy.  A simple  procedure  for  elimination  of  some  of  these 
constants  is  to  estimate  the  real  parts  of  the  frequenices  (for 
1=0)  dominated  by  these  modes  and  compare  them  with  the  highest 
flexural  frequency.  Based  on  this  comparison,  the  relative  im- 
portance of  these  terms  can  be  easily  determined.  The  rows  and 
columns  corresponding  to  the  terms  which  are  not  important  may 
then  be  omitted  from  the  system  of  equations  in  (C5  ). 

Panel  Flutter  of  [0/90]  Laminate 


For  a laminate  made  of  orthotropic  layers  with  their  material 

axes  parallel  to  the  x and  y directions,  the  solutions  associated 

with  V®  and  U/  are  uncoupled.  Hence,  the  constants  d and  e„ 
y n n 

(n  = 1,2...N)  can  be  removed  from  equation  (C 5 ) . Also,  b^  and 
(n  = 1,2...N)  can  be  eliminated  from  equation  (C 5 ) provided 
that 
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(A,  , a 
11  n 


— 2 2 ^ 
P-  ) <°ll«n  ^ ^55 


- Iw^)  - 
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R(d^)^  ^ 0, 
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This  process  of  elimination  yields 


N } + iSiD  - Poj^la 

A XX  ^ n 

n 


a a 2 

+ 41  E ^-iL- 

a - a 

m=l  n m 
n-m=odd 


= 0 , n = 1,2, ..  .N 


(Cll) 


Equation  (Cll)  is  a system  of  N homogeneous  equations  and  for  a 
nontrivial  solution,  the  determinant  of  the  coefficients  a^  must 
vanish.  The  critical  flutter  parameter,  may  be  determined 

as  the  lowest  value  of  1 for  which  the  real  and  imaginary  parts 
of  the  determinant  are  zero.  The  determinant  is  evaluated  for  a 
range  of  values  of  real  w and  A to  determine  these  quantities 
at  which  both  the  real  and  imaginary  parts  of  the  determinant 
vanish. 
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APPENDIX  D. 

SOLUTION  OF  THE  MOISTURE  DIFFUSION  PROBLEM 


The  governing  differential  equation  of  the  one-dimensional 
moisture  diffusion  equation  is 


at  - -53^  -5^) 


(Dl) 


where  C is  the  moisture  concentration.  If  the  derivatives  with 
respect  to  space  coordinate  x are  replaced  by  finite  differences, 
the  partial  differential  equation  (Dl)  can  be  converted  to  a 
set  of  ordinary  differential  equations 


ac. 


Tt-  - ^i-l^i-1  + ‘^i'^i  + ‘^i+l^i+i 


1 — 2,3,...  n 


(D2) 


where 


«i_l  = (4Mi  - + yi_i)/4d' 

,2 


6^  = -2Mi/d' 


{D3) 


i.  = (4ui  + - y._;^)/4d' 


= diffusivity  at  node  i 


d = mesh  size  = H/n  when  all  meshes  are  of  equal  size, 
H is  the  laminate  thickness  and  n+1  = total  number  of  nodes 
including  the  boundary  points. 


With  prescribed  boundary  conditions  and  as  well  as 


initial  distribution  of  C known,  the  set  of  equations  (D2)  can 
be  easily  solved  by  a predictor  corrector  scheme.  With  the 
use  of  this  scheme,  the  dif fusivities  may  be  considered  to 


be  dependent  on  C^, 
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To  consider  jump  discontinuities  in  dif fusivities,  the 
equation  (D2)  at  the  node  k which  coincides  with  the  interface 
between  layers  m and  m+1,  is  modified  by  the  use  of  the  following 
values  of  *^i+l' 


6 = (y"'  + yj  )/[d  (d^  + d )] 

i~i  k k— 1 m m m+1 


«i  = -[(<  + + ll-r"  + 


where  d^^  and  are  the  mesh  sizes  in  layers  m and  m+1,  respec- 

tively. These  values  are  obtained  on  the  assumption  that  the 
moisture  concentration  at  node  k is  the  same  in  layers  m and 
m+1.  The  following  relations  have  been  made  use  of  to  obtain 
equation  (D4) : 


^m  TT  ^ ^^k  ■ ^k-1^  ^^k  ~ ^k-l^^'^m  ^k  ^^k  " ^S-1 


‘^m+l  TT  = (C^+1  - Ck)/d^+l  + C 

^^k  ■ ^^k+1  ^ ^k+2^/‘^m+l 

“k  - "S-l  * =k-2>/2<3„  = C' 


■^k+2^/2'^m+l 


Equations  (D5)  and  (D6)  are  obtained  from  (Dl)  by  the  use 
of  backward  and  forward  difference  formulae  at  node  k in  layers 
m and  m+1,  respectively.  A mass  conservation  condition  (mois- 

m ^*"k  m+1 

ture  transfer)  across  the  interface,  i.e.  u,  — .r — = M,  — 

k 9x  k 9x 

in  the  form  of  a higher  order  difference  formula  yields  the  rela- 
tion (D7) . The  addition  of  equations  (D5)  and  (D6)  and  the 
use  of  (D7)  result  in  (D4). 

In  the  present  form,  the  computer  program  solves  the  system 
of  equations  (D2)  - (D4)  for  prescribed  initial  and  boundary  con- 
ditions as  well  as  a given  temperature  distribution.  Diffusivi- 
ties  at  various  temperatures  and  moisture  concentrations  are 
input  in  a tubular  form  and  linear  interpolation  is  used  to  cal- 
culate diffusivity  for  any  temperature  and  moisture  concentra- 
tion. 
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Epoxy  Matrix 


Table  2.  Various  Damping  Levels  for  Epoxy  Matrix 


Level  of  Damping 


2 (Base;  table  1) 


Multiplication  Factor 
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0.25 


2.5 


r 


Table  3.  Effect  of  Damping  Levels  on  Flutter  Parameter 

X for  Laminates  1,  2,  and  3. 
cr 


Laminate 

JL. 

Boron/epoxy , 

[0/90]g,  a = 

0.6  in , 

H = 0.06  in 

Level  of 
Damping 

Figure  No. 
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in 

Level  of 
Damping 

0 

Figure 
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cr 

30 
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31 

130 
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Table  4. 


Laminate 


Laminate 


Level 


Effect  of  Damping  Levels  on  Flutter  Parameter 
for  Laminates  4 and  5. 


4.  Graphite/Epoxy  [0/9012^/  a = 3.6  in.,  H = .06  in 


Level  of  Damping  ^cr 

0 340 

1 157 

2 150 

3 150 
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Imaginary 


Figure  3.  Variation  of  Lamina  and  Matrix  Complex  Young's  Moduli  With  Frequency 
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Figure  10.  Shear  Correction  I 
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Figure  11.  Shear  Correction  Factors 
f±15°l  Laminate 


and  k 
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Figure  13.  Shear  Correction  Factors  ^44  and  for 

[±30®lg  Laminate 
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Figure  14. 
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Shear  Correction  Factor  k'*  for  t±30®]  Laminate 
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Figure  21.  Variation  of  Shear  Correction  Factors  With  the 

Fraction  of  ±45°  Layers  in  a [0/±45/90]g  Laminate 
(Equal  Fractions  of  0°  and  90°  Layers) 
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Figure  22.  Variation  of  Shear  Correction  Factors  With  the 
Fraction  of  ±45°  Layers  in  a [0/±45]  Laminate 
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Figure  24.  Variation  of  Shear  Correction  Factors  With  the 
Fraction  of  ±45®  Layers  in  a f±45/0]  Laminate 
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9 , Degrees 

Figure  25.  Shear  Correction  Factors  k^^,  k^^  for  Antisymmetric 
[±9]  Laminates 


COATING 


Figure  26.  Moisture  Distribution  in  a Coated  Composite 
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Figure  31.  Flutter  Parameter  - Frequency  Plot 


220 


Figure  32.  Flutter  Parameter  - Frequency  Plot 
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Figure  33.  Flutter  Parameter  - Frequency  Plot 


Figure  34 . Flutter  Parameter  - Frequency  Plot 
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Figure  36.  Flutter  Parameter  - Frequency  Plot 


Figure  37.  Flutter  Parameter  - Frequency  Plot 
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Figure  40.  Frequency  Spectrum  for  [0/90/0/90]  Cross  Ply  Laminate  5,  a/h 


Effect  of  Environment,  Aerodynamic  Damping  and  Prestress  on  Amplitude 
Ratio  VS-  > for  (±45/0, /90]  _ Laminate  6,  a/h  = 40,  h = .06  in. 


